Conpackt, A Contouring Package for Triangular Meshes

David J. Kennison
NCAR, P.O. Box 3000, Boulder, Colorado 80307-3000
email: kennison@ucar.edu

Table of Contents

INTRODUCTION

This document describes an NCAR Graphics package called CONPACKT that allows a user to construct contour plots from data on triangular meshes. CONPACKT provides a sort of tool kit of FORTRAN subroutines that can be called in various combinations to draw different kinds of contour plots.

This section is intended to give an overall view of CONPACKT and selected aspects of its design; it covers some details, but, in general, one should refer to the sections "SUBROUTINES" and "PARAMETERS" for detailed descriptions of subroutines and parameters mentioned. (Parameters are mentioned by name; all the names are of the form 'XXX', where XXX is a three-character mnemonic.) The section "ERROR HANDLING" describes error messages written by CONPACKT. The section "EXAMPLES" describes the examples available for CONPACKT.

It is assumed that the reader is familiar with NCAR Graphics in general and has some knowledge of the packages AREAS and EZMAP/EZMAPA.

Triangular Mesh Structure

Note: See the example "ctex01" (both the code and the first frame that it produces) for help in understanding the structure of a triangular mesh.

To represent a triangular mesh requires three singly-dimensioned arrays: RPNT defines points, IEDG defines edges (each of which consists of two connected points), and ITRI defines triangles (each of which consists of three connected edges). The elements of each array form "nodes" having nominal lengths as follows:

      PARAMETER (LOPN=4)  !  length of a point node
      PARAMETER (LOEN=5)  !  length of an edge node
      PARAMETER (LOTN=4)  !  length of a triangle node
(In some cases, additional elements are needed in some types of nodes, but these are the basic ones that must be present.)

The four elements of a point node are reals, as follows:

  1. the X coordinate of the point;
  2. the Y coordinate of the point;
  3. the Z coordinate of the point;
  4. the field value at the point;
The five elements of an edge node are integers, as follows:

  1. the base index, in RPNT, of the node defining point 1 of the edge;
  2. the base index, in RPNT, of the node defining point 2 of the edge;
  3. the base index, in ITRI, of the node defining the triangle to the left of the edge, plus the index (1, 2, or 3) of the edge within that triangle node (-1 if there is no triangle to the left);
  4. the base index, in ITRI, of the node defining the triangle to the right of the edge, plus the index (1, 2, or 3) of the edge within that triangle node (-1 if there is no triangle to the right);
  5. a utility flag for use by algorithms that scan the structure.
The four elements of a triangle node are integers, as follows:

  1. the base index, in IEDG, of the node defining edge 1 of the triangle;
  2. the base index, in IEDG, of the node defining edge 2 of the triangle;
  3. the base index, in IEDG, of the node defining edge 3 of the triangle;
  4. a flag set non-zero to block use of the triangle, effectively removing it from the mesh. For simple purposes, the values 0 and 1 suffice; however, see the routines CTTDBF and CTTDBM and the internal parameters 'TBX' and 'TBA' for a description of a scheme allowing one to selectively block triangles by using different bits of this flag.
The "left" and "right" sides of an edge are defined as they would be by an observer standing on the front side of the surface defined by the triangular mesh, at point 1, and looking toward point 2, of the edge. It is possible, if there are "holes" in the mesh, that there will be no triangle on the left or on the right of an edge, but there must be a triangle on one side or the other.

Note: The edges pointed to by a particular triangle node must be given in counter-clockwise order, as viewed from the chosen "front" side of the mesh. In fact, it is the ordering of the nodes that defines which side is the "front" and which side is the "back". In the case of a sphere, one would probably use the outside of the sphere as the "front". In the case of a Moebius strip (which I have experimented with a bit), there is a problem with this - you either have to have a seam across the strip or go around it twice to avoid having a seam - but that's probably not a case of great practical interest ... :-)

The "base index" of a point node, an edge node, or a triangle node is always a non-negative multiple of the length of the node, to which can be added an offset to get the index of a particular element of the node. For example, if IPTT is the base index of a triangle of interest, ITRI(IPTT+1) is its first element, which is the base index of the triangle's first edge. Thus, IEDG(ITRI(IPTT+1)+2) is the second element of the triangle's first edge; that is to say, it is the base index of the second point of the first edge of the triangle with base index IPTT. In a similar fashion, it may be seen that RPNT(IEDG(ITRI(IPTT+1)+2)+3) is the third (Z) coordinate of the second point of the first edge of the triangle with base index IPTT.

It is the pointers from the edge nodes back to the triangle nodes that allow CONPACKT to navigate the mesh, moving from triangle to triangle as it follows a contour line. These pointers are a little tricky to define: if IPTE is the base index of an edge node and IEDG(IPTE+3) is zero or more, saying that there is a triangle to the left of the edge, then IEDG(IPTE+3) is the actual index of that element of the triangle node that points to the edge node; that is, ITRI(IEDG(IPTE+3))=IPTE. The base index of the triangle node defining that triangle is IPTT, where IPTT=LOTN*((IEDG(IPTE+3)-1)/LOTN), and the index of the pointer to the edge within the triangle node is IPTI=IEDG(IPTE+3)-IPTT, so that ITRI(IPTT+IPTI)=IPTE. Similar comments apply to element 4 of an edge node, which points into the triangle node defining the triangle to the right of the edge.

For some types of triangular mesh, the maximum number of points, edges, and triangles can be computed easily:

Once we know, at most, how many points, edges, and triangles we're going to have, we can set parameters defining the space to be reserved for the triangular mesh:

      PARAMETER (MPNT=MNOP*LOPN)  !  space for points
      PARAMETER (MEDG=MNOE*LOEN)  !  space for edges
      PARAMETER (MTRI=MNOT*LOTN)  !  space for triangles
Then, we can declare the arrays to hold the point nodes, edge nodes, and triangle nodes defining the triangular mesh:

      DIMENSION RPNT(MPNT),IEDG(MEDG),ITRI(MTRI)
In all of the examples, code like that above (throughout the section) will be found.

Types of Triangular Meshes

The triangular meshes that I have dealt with and provide examples for might be classified as follows:

In the plane:

"ctex01" and "ctterr": These are the examples to look at if you want to learn the details of constructing a triangular mesh from a rectangular grid in the plane completely from scratch. "ctex01" uses a small grid (7 x 7, by default), while "ctterr" is an adapted version of the code that uses a rectangular grid contributed by a user, representing a vertical slice through an ocean basin, and deformed to follow the terrain of the ocean bottom; this grid is dimensioned 4000 x 62, so the triangular mesh consists of 3999 x 61 x 2 = 487,878 triangles, which is a pretty good workout for CONPACKT. Both examples contain a nice little routine, called DRWMSH, to draw the mesh, but, in the interest of making the metafile a bit smaller, the call to it in "ctterr" has been commented out. "ctex01" also contains code to label the points, edges, and triangles of the mesh, but these have been omitted from "ctterr".

"ctex02": Demonstrates a possible way to construct, in the plane, a triangular mesh that fills a particular area of interest (in this case, the contiguous "lower 48" states of the US). Again, this example constructs the triangular mesh completely from scratch; it begins with a set of equilateral triangles that completely fill the map space, discards those that are wholly outside the area of interest, and adapts the others to better fit that area.

In three dimensions, on the globe, rectangular:

Each of these examples uses EZMAP routines to show the earth in orthographic projection, as seen from four different viewpoints, and each calls the routine CTTMRG to deal with a rectangular grid mapped onto the surface of the globe. In the code for any of these examples, search for the call to CTTMRG to see how the triangular mesh is constructed.

"ctllg1": In this example, the grid is one created for illustrative purposes; it is dimensioned 81 x 41 and is mapped onto the portion of the globe between longitudes 80W and 80E and latitudes 40S and 40N, so that it covers a limited portion of the globe.

"ctllg2": In this example, the grid is one created for illustrative purposes; it is dimensioned 171 x 81 and is mapped onto the portion of the globe between longitudes 170W and 170E and latitudes 80S and 80N, so that it nearly covers the globe.

"ctllg3": In this example, the grid is one created for illustrative purposes; it is dimensioned 181 x 91 and is mapped onto the entire globe. All of the points along the top edge of the rectangular grid map to the North Pole and all the points along its bottom edge map to the South Pole; the left and right edges of the rectangular grid meet along a seam that runs from the North Pole to the South Pole along the International Date Line.

"ctgaus": This example uses a Gaussian grid from Paul Swarztrauber; it is dimensioned 129 x 257 and is rotated 90 degrees as compared to the grids of examples "ctllg1", "ctllg2", and "ctllg3". It is mapped onto most of the globe. The left and right edges of the grid map to little circles around the North and South poles, respectively, while the top and bottom edges of the grid meet along a seam that runs along a meridian from pole to pole.

"ctpopg": This example uses a typical POP grid from a user site; it is dimensioned 101 x 116 and is mapped onto the entire globe, less two little circles. Think of this grid as being wrapped around the globe so that its left and right edges meet along a seam on a meridian and then elastically deformed so as to pull the little circle formed by its top edge to a position over Greenland, while leaving the little circle formed by its bottom edge centered over the South Pole and wholly contained within Antarctica. Finally, portions of the grid over land are omitted, as this grid is used for ocean modelling.

"ctorca": This example uses a typical ORCA grid from a user site; it is dimensioned 181 x 148 and is mapped onto the entire globe, less a circle at the South Pole and a few other areas over Asia and Africa. Think of wrapping this grid around the globe so that its left and right edges meet along a seam on a meridian and then squashing flat the circle formed by its top edge in such a way that it becomes an extension of the left/right seam passing through the North Pole. (The circle formed by the bottom edge of the grid is centered on the South Pole and remains open.) Certain rectangular portions of the grid are then severed from the rest of it and shifted elastically so as to provide denser coverage of certain areas (the Mediterranean, Black, Caspian, and Red Seas). Finally, portions of the grid over land are omitted, as this grid is used for ocean modelling.

"ctfite": This example uses a POP grid from a user (I got it from Fred Clare); it is dimensioned 321 x 384 and is mapped onto the globe in much the same way as the grid of "ctpopg". The data are actual user data and show ice thickness, so there is little of interest except at the poles.

"ctswth": This example uses a rectangular grid from Simona Bordoni and represents the area seen by a satellite during one orbit of the earth. I think the data are wind speeds over the ocean.

In three dimensions, on the globe, triangular:

Each of these examples uses EZMAP routines to show a portion of the earth in orthographic projection. Each calls the routine CTTMTL to deal with an intrinsically triangular mesh mapped onto the surface of the globe. In the code of any of these examples, search for the call to CTTMTL to see how a triangular mesh of the form expected by CONPACKT is constructed from triangles provided by a user.

"ctcbay": This example uses data provided by Tom Gross. It consists of the latitudes and longitudes of 7,258 points on the globe, plus an array of 13,044 index triplets, each specifying the indices of three of those points forming a triangle. The triangular mesh represents Chesapeake Bay.

"ctnccl": This example uses data provided by Brett Estrade. It consists of the latitudes and longitudes of 32,218 points on the globe, plus an array of 58,641 index triplets, each specifying the indices of three of those points forming a triangle. The triangular mesh represents the ocean off the coast of North Carolina.

In three dimensions, on the globe, geodesic:

Each of the examples "ctgeo1", "ctgeo2", and "ctgeo3" uses EZMAP routines to show the earth in orthographic projection, as seen from four different viewpoints, and each uses a geodesic grid constructed using a different method. In all three cases, the net effect is to subdivide the faces of an icosahedron (a regular Platonic solid with twenty equilateral triangles as faces), into smaller triangles and then project the vertices of those triangles outward onto the surface of a circumscribed sphere. Subdivision of a triangle is done using N-1 lines parallel (in some cases, roughly parallel, I think) to each of its three sides, dividing each side into N pieces; by default, N = 16, but, in "ctgeo1" and "ctgeo2", this is set by a PARAMETER statement and so is easily changed. Each face is thereby subdivided into N**2 (by default, 256) smaller triangles, so the total number of triangles is 20*N**2 (by default, 5120). The number of edges is 30*N**2 (by default, 7680) and the number of vertices is 10*N**2+2 (by default, 2562).

"ctgeo1": In this example, the geodesic grid is constructed in a clever way that makes very efficient use of memory, but the code to do it, in the subroutine GTGEO1, will probably prove quite difficult to follow. Each triangular face of the icosahedron is subdivided, in a single step, into N**2 smaller triangles of the same size. Once those triangles are projected outward onto the sphere, some of them (those near the center of a face of the icosahedron) are quite a bit larger than others (those near a vertex of the icosahedron).

"ctgeo2": In this example, two innovations were introduced in the subroutine GTGEO2, which constructs the geodesic grid:

"ctgeo3": This example makes use of data from the NetCDF file "C02562.global.nc", which was downloaded from the following Web site:

Go to "Step 2", click on "sample spherical geodesic grids", and then on "Geodesic grid with 2562 cells". (Actually, "ctgeo3" can easily be modified to use either of the other grids mentioned there, as well.)

These data are of interest in two ways:

"ctgc23": This example produces four plots comparing the geodesic grids of examples "ctgeo2" and "ctgeo3". Two of the plots show the grids superimposed (on a hemisphere of the globe and in a relatively close-up view); the other two plots show histograms of the areas of the triangles produced by the two methods.

In three dimensions, on the globe, other:

"ctiscp": This example uses a grid created by the "International Satellite Cloud Climatology Project" (ISCCP). It consists of 6596 cells on the globe, each of which is a rectangle on a cylindrical equidistant projection. The cells are created by dividing the surface of the earth into strips, parallel to the equator, each of which is 2.5 degrees of latitude wide, and then subdividing each strip into an integral number of cells in such a way that all of the cells have almost the same surface area on the globe. The subroutine GTISCP triangulates this grid by connecting the centers of "adjacent" cells and then using the connecting lines to form triangles. The code is rather difficult to describe and may be a bit difficult to read.

"ctisc2": This example represents an attempt to deal with the ISCCP grid in a more understandable way. It treats each cell as a polygonal patch like those of example "ctgeo3" and then uses the technique of that example to determine which of the polygonal patches are adjacent. (Some of the polygonal patches have one or more sides of zero length.) Whether the resulting code is easier to understand than that of "ctiscp" is debatable.

"ctwng1" and "ctwng2": These examples use SEAM grids (where SEAM stands for "Spectral Element Atmospheric Model") obtained from Houjun Wang. See the following Web site for more informataion about such grids:

The basic elements of each grid are quadrilateral patches on the surface of the globe, obtained by subdividing the faces of a cube and projecting the pieces outward onto the surface of the globe. The patches are further subdivided into smaller quadrilaterals; in "ctwng2", in particular, this is done in such a way as to provide higher-resolution coverage of the US. Each of the quadrilaterals is then split in half along a diagonal to form two triangles and the routine CTTMTL or CTTMTX is used to form a triangular mesh of the form expected by CONPACKT routines.

In three dimensions, on an arbitrary surface:

"cttd01" and "cttd02": These examples demonstrate that CONPACKT is not limited to dealing with triangular meshes in a plane or on the surface of the globe. In each case, the triangular mesh used represents an arbitrary surface in 3-space (looking, perhaps, a bit like one of the dancing mushrooms from Disney's "Fantasia"); one of the meshes forms a closed object, while the other has a hole in it through which the "back" side of the mesh is visible. The 3D package TDPACK is used to draw perspective views of these grids, with contours on them.

The routine CTTDCA is used to create a cell array that creates an accurate "picture" of the object. The routine CTTDBF is used to set bits in the blocking flags for the triangles to reflect various conditions (triangle hidden behind another, triangle seen from the back side, triangle seen nearly edge-on) and CTTDBM is used, just before drawing the mesh itself or just before drawing contours on the mesh, to allow drawing in only those triangles that, for example, are seen from the front, are not hidden, and are not too nearly edge on. The net result is to solve the "hidden-line" problem for the mesh or for the contours on a triangle-by-triangle basis (not quite as good as one would like, but not too bad, either).

Routines Which May Be Called

The following routines are meant to facilitate the process of creating a triangular mesh of data:

Once a triangular mesh has been defined, the code to draw a contour plot will probably begin with several calls to set internal parameters affecting the behavior of the routines called after that. All the internal parameters have default values; only those which are to have values different from the default need to be set. Routines which may be called to set the values of parameters are as follows:

In general, once a parameter is given a value by a call to one of these routines, it retains that value until a similar call resets it. Thus, many of the parameter-setting calls need to be done only once and do not need to be repeated for each new contour plot.

After all required parameters have been set, the process of drawing a contour plot begins with a call to an initialization routine:

This routine CTMESH performs initialization required for subsequent calls to work properly. The dimensions of the arrays defining the triangular mesh and the dimensions of the real and integer workspace arrays are copied to internal variables. Pointers in the array defining the edges of the mesh are adjusted so as to make the data value at the first point of each edge less than or equal to the data value at the second point of each edge, which is important for various algorithms in other routines. Pointers used to manage the real and integer workspaces are initialized. It is decided how the 3D coordinates of the points of the triangular mesh are to be mapped into 2D user coordinates and how those coordinates are to be mapped to the plotter frame. If the user has not called the SPPS routine SET (or done the equivalent GKS calls), that call is done. The minimum and maximum values in the data array are found and it is decided whether or not the data field is essentially constant. If contour levels are to be chosen automatically, the parameter arrays specifying the contour levels are cleared, so that automatic selection will take place when these levels are required.

Among the internal parameters are arrays completely specifying contour levels of interest and what is to be done at each of those levels. These arrays may be used to take complete control of the contouring process; most users will probably elect not to do this, but to let CONPACKT choose the levels. At this point, then, as a rule, none of these parameter arrays will have been filled. No calls need be done to fill them; if they are empty when they are needed, the required values will be chosen at that point. For certain applications, however, it is desirable to force the selection of contour levels and perhaps the character strings which are to be used as contour labels. This may be done by means of calls to one or both of the following routines:

The advantage of calling one or both of these routines is that, after each call, one can modify what the routine called has done in order to produce desired effects. For example, after CTPKCL is called, one might set elements in parameter arrays which will cause the labelled contour lines to be solid and the unlabelled contour lines to be dashed. As another example, after CTPKLB is called, one might check for the string "0" as a label and change it to "0.0".

At this point, various other routines may be called:

Four routines may be called when 'MAP' is set to 2, which says that the 3D package TDPACK is being used to project the triangular mesh into the plane.

Finally, to advance the frame, the user must call the SPPS routine FRAME; CONPACKT won't do it.

At any time, it is possible to retrieve the value of an internal parameter by calling one of the three following routines:

Several routines in CONPACKT are not called by the user, but by CONPACKT itself. The default versions of these routines, in all cases but one, do nothing; the routines exist simply to be replaced by the user. These routines are as follows:

Two additional routines are provided for use in error-recovery situations:

Different Styles of Contour Plots

The design of CONPACKT allows one to construct contour plots in many different styles. Which style is chosen will depend on the amount of computer time and memory available, on the capabilities of the device on which the plots are being drawn, and on the intended use of the plots. Some possibilities are discussed below.

If the intent is to draw the contour plot using relatively few computer resources, the following sequence of calls will suffice:

  1. Call parameter-setting routines.

  2. Call CTMESH to initialize the drawing of the contour plot.

  3. Call CTBACK to draw a simple background.

  4. Call CTCLDR to draw the contour lines, with labels generated by a dashed-line package.

  5. Call CTLBDR to draw an informational label below the plot and high/low labels on the plot.

Assuming that step 1 is null and that all default parameter values are used, the resulting plot will look much like what would have been drawn by CONREC. The labels won't be positioned very well and lines may be drawn through them. The characters used for the informational label and for high and low labels will be of better quality than those used for the line labels; if desired, additional resources may be saved by a parameter-setting call forcing the PLOTCHAR package to use lower-quality characters.

To produce a better-positioned set of labels, step 1 above should include a call to CTSETI setting the parameter 'LLP', which says how line labels are to be positioned, to one of the values 3 or -3; this turns off the generation of line labels by a dashed-line package called by CTCLDR and causes them to be positioned and written by the routine CTLBDR instead, using a penalty scheme which produces quite pleasing results. Contour lines will still pass through line labels positioned in this way.

Contour lines may be prevented from passing through labels drawn by CTLBDR in one of two ways, depending on the nature of the plotting device. If the device allows for the use of the GKS fill area primitive and if the result of drawing one object on top of another is that the pixels affected by drawing the second object simply change color (as happens on most terminals, but not usually on a device which exposes film, for example), then one may insert calls changing the values of the parameters 'ILB', 'HLB', and 'LLB' in such a way as to force boxes surrounding the informational label, the high and low labels, and the line labels to be filled with the background color prior to the drawing of the labels. This has the effect of preventing the contour lines from passing through the labels.

If filling the label boxes will not work, then a software technique, using routines in the utility AREAS, may be used instead. The sequence of calls will then be as follows:

  1. Call parameter-setting routines.

  2. Call CTMESH to initialize the drawing of the contour plot.

  3. Call CTBACK to draw a simple background.

  4. Call ARINAM (in AREAS) to initialize an area map.

  5. Call CTLBAM to generate a list of label positions and to put label boxes for the labels at those positions into the area map.

  6. Call CTCLDM to draw contour lines masked against the area map. Each line is broken into pieces lying outside label boxes and pieces lying inside label boxes and only the former are drawn.

  7. Call CTLBDR to draw the labels in the boxes left vacant by step 6.

Some users may wish to draw a solid-filled contour plot, rather than a line plot. This is also done using routines in the utility AREAS, as follows:

  1. Call parameter-setting routines.

  2. Call CTMESH to initialize the drawing of the contour plot.

  3. Call ARINAM (in AREAS) to initialize an area map.

  4. Call CTCLAM to put contour lines into the area map.

  5. Call ARSCAM (in AREAS) to scan the area map and to recover from it the polygons created by the lines in the area map. A user-defined routine is called to fill (or not to fill) each polygon. Filling may be done by calls to routines in the utility SOFTFILL or to the GKS routine GFA.

Labels may be written on a solid-filled contour plot; this requires adding a call to CTLBAM after step 4 and a call to CTLBDR after step 5. It may be better, though, to draw a key for the solid-filled plot. There is no routine in CONPACKT to do this, but the NCAR Graphics utility LABELBAR provides such a capability.

The packages AREAS, CONPACKT, and EZMAP/EZMAPA may be used cooperatively to achieve other desired effects. For example, if the contour plot being drawn represents output from an ocean model, it may be desirable to draw contours (or to fill contour bands) only over the oceans on a background drawn by EZMAP.

A solid-filled contour plot may be drawn in another way, using the GKS primitive "cell array". The steps required are as follows:

  1. Call parameter-setting routines.

  2. Call CTMESH to initialize the drawing of the contour plot.

  3. Call CTCICA to generate color indices in a cell array.

  4. Call the GKS routine GCA to "paint" the cell array.

Using cell arrays has advantages and disadvantages compared to using AREAS to generate GKS "fill areas". One disadvantage is that, for the same sort of quality, a rather high-resolution cell array must be used; it may be time-consuming and space-consuming to generate and to display such a cell array; in particular, using the "zoom" capability of "idt" on a cell array doesn't work too well. Another disadvantage of the cell-array approach is that it is more difficult (sometimes impossible) to limit the coloring of contour bands to specific areas (as, for example, over land or over ocean), which can be done relatively easily using the fill-area approach. The major advantage of the cell-array approach is that it avoids the problems that, to some extent, accompany the use of the package AREAS. (While many of these problems have been fixed, one can still expect an occasional glitch, particularly when using certain EZMAP projections.)

Contour Level Selection

The routines CTCLAM, CTCLDM, CTCLDR, CTLBAM, and CTLBDR all generate output for specified sets of contour levels. The parameter 'NCL' specifies the number of different contour levels for which something is to be done by one or more of the above routines. The value of 'NCL' may not exceed 256. For each value of I from 1 through 'NCL', the Ith element of the parameter array 'CLV' specifies a contour level and the Ith element of each of the parameter arrays 'AIA', 'AIB', 'CLC', 'CLD', 'CLL', 'CLU', 'LLC', and 'LLT' is an associated datum, saying something about what is to be done at that contour level. Detailed descriptions of these parameter arrays are given in the section "PARAMETERS", below; thumbnail descriptions of elements of these parameter arrays are as follows:

The contents of these parameter arrays may be specified completely by CONPACKT, completely by the user, or some of both. The parameter 'CLS' says whether or not CONPACKT is to select contour levels and associated quantities and, if so, in what manner. The default value of 'CLS' indicates that CONPACKT is to select about 16 contour levels at multiples of some nice interval, which it is to choose. By default, then, the call to CTMESH zeroes 'NCL'; during the first subsequent call which requires contour levels to have been chosen, they are chosen. The associated parameter elements are set to indicate which contour lines should be labelled, what character strings should be used as labels, which lines should be put in an area map by a call to CTCLAM, and what area identifiers should be used for areas "above" and "below" these lines.

Other values of 'CLS' may be used to change the way in which CONPACKT chooses contour levels. See the detailed descriptions of 'CLS' and of the other parameters 'CIS', 'CIT', 'CIU', 'CMN', 'CMX', 'LIS', 'LIT', and 'LIU'.

The user may elect to set all of the contour levels and associated quantities. Suppose, for example, that it is desired to draw labelled solid lines for each of the values .1, .2, .3, ..., .9 and unlabelled dashed lines for each of the values .05, .15, .25, ... .95. The following code, inserted before the call to CTMESH, will set the required parameters:

      CALL CTSETI ('CLS - CONTOUR LEVEL SELECTION',0)
      CALL CTSETI ('NCL - NUMBER OF CONTOUR LEVELS',19)
      DO 101 I=1,19
        CALL CTSETI ('PAI - PARAMETER ARRAY INDEX',I)
        CALL CTSETR ('CLV - CONTOUR LEVEL VALUE',REAL(I)/20.)
        IF (MOD(I,2).EQ.1) THEN
          CALL CTSETI ('CLU - CONTOUR LEVEL USE',1)
          CALL CTSETI ('CLD - CONTOUR LINE DASH PATTERN',21845)
        ELSE
          CALL CTSETI ('CLU - CONTOUR LEVEL USE',3)
          CALL CTSETI ('CLD - CONTOUR LINE DASH PATTERN',65535)
        END IF
  101 CONTINUE
In the above code, 'CLS' is zeroed to suspend the selection of contour levels by CONPACKT itself. Then, 'NCL' is set to say how many contour levels are to be defined. Then, in a loop on I from 1 to 19, 'PAI' is set to tell CONPACKT which element of each parameter array is to be set, the Ith element of 'CLV' is set to REAL(I)/20., which, for each I, gives one of the desired contour levels, the Ith element of 'CLU' is set to a 1 if just the line is to be drawn or to a 3 if both the line and the labels for the line are to be drawn, and the Ith element of 'CLD' is set to 21845 (octal 52525) if a dashed line is to be used or to 65535 (octal 177777) if a solid line is to be used.

Note that 'NCL' must be set prior to setting any element of 'CLV' or the associated arrays.

Note also that, when an element of 'CLV' is set, all of the associated elements of the associated arrays receive a default value. (In fact,the default element of 'CLU' is 1, and the default element of 'CLD' is a pattern specifying a solid line, so two of the calls in the code above are redundant.)

Note also that the use of CTSETI to specify the value of 'CLD' will work only if 'DPU' is positive, implying the use of DASHCHAR routines; if 'DPU' is negative, DASHPACK routines will be used and, by default, they will expect underscores, rather than apostrophes, to represent gaps, so that CTSETC calls will be needed to define the dash patterns.

In some cases, the user will want to let CONPACKT choose a set of contour levels and then either add other levels of interest, modify elements of the associated parameter arrays, or both. Suppose, for example, that it is desired to have CONPACKT pick the levels, that contour lines at positive levels are to be drawn in red, that contour lines at negative levels are to be drawn in blue, and that contour lines at the zero level are to be drawn in white. The following code, inserted after the call to CTMESH, would do the job (assuming that color indices IBLU, IRED, and IWHI have previously been defined):

      CALL CTPKCL (...)
      CALL CTGETI ('NCL - NUMBER OF CONTOUR LINES',NOCL)
      DO 101 I=1,NOCL
        CALL CTSETI ('PAI - PARAMETER ARRAY INDEX',I)
        CALL CTGETR ('CLV - CONTOUR LEVEL VALUE',CLEV)
        IF (CLEV.LT.0.) THEN
          CALL CTSETI ('CLC - CONTOUR LINE COLOR INDEX',IBLU)
        ELSE IF (CLEV.GT.0.) THEN
          CALL CTSETI ('CLC - CONTOUR LINE COLOR INDEX',IRED)
        ELSE
          CALL CTSETI ('CLC - CONTOUR LINE COLOR INDEX',IWHI)
        END IF
  101 CONTINUE
In the code above, the routine CTPKCL is called to force CONPACKT to pick a set of contour levels. Then, the value of 'NCL' that it chose is retrieved, and a loop is run from 1 to that value. On each pass through the loop, the parameter array index 'PAI' is set to tell CONPACKT what element of the parameter arrays is being accessed and then one of the contour levels chosen is retrieved to the variable CLEV. Depending on the value of CLEV, the associated element of the parameter array which specifies the color index of the contour lines at that level is set to produce a blue line, a red line, or a white line.

CONPACKT Coordinate Systems

The mapping of contours into the plotter frame depends on three things:

By default, to describe contour lines and other objects on the contour plot, CONPACKT generates X, Y, and Z coordinates representing points on the edges of the triangles defined by the data of the triangular mesh and then discards the Z coordinate, using only the X and Y coordinates as the user coordinates in calls to NCAR Graphics routines.

If the internal parameter 'MAP' is given a positive non-zero value, each triplet of X, Y, and Z coordinates is mapped, prior to use, by a statement of the form

      CALL CTMXYZ (IMAP,XINP,YINP,ZINP,XOTP,YOTP)
IMAP is the value of 'MAP'; XINP, YINP, and ZINP are the unmapped (input) 3D coordinates; and XOTP and YOTP are the mapped (output) 2D coordinates to be used as "user coordinates".

The default version of CTMXYZ does the following mappings (where RTOD = 57.2957795130823 (180/pi):

The subroutine CTMXYZ may be replaced by a version which does other desired mappings. If the CONPACKT routines are loaded from a binary library, this can usually be done by just compiling one's own version of the routine, so that it replaces the one from the library. The definition of mapping 1 should not be changed.

The way in which "user coordinates" are mapped to "fractional coordinates" in the plotter frame is determined by the current definitions of the "window" in the user system and the "viewport" on the plotter frame. The window and viewport may have been defined by a call to the SPPS routine SET or by calls to GKS routines; the former will be described.

A call to the SPPS routine SET has the form

      CALL SET (XVPL,XVPR,YVPB,YVPT,XWDL,XWDR,YWDB,YWDT,LNLG)
All arguments are REALs except for LNLG, which is an INTEGER. The first four arguments must all be between 0 and 1, inclusive; they define a rectangular area in the fractional coordinate space of the plotter frame known as the "viewport". The next four arguments define a rectangular area in "user coordinate" space known as the "window". The final argument indicates whether the mapping of user coordinates into the viewport is to be linear or logarithmic in X and Y. See the documentation of the package SPPS for further details.

By default, CONPACKT (specifically, the routine CTMESH) calls SET. One may, by setting the parameter 'SET' to zero, prevent CONPACKT from doing this; in that case, one must do the call for oneself or depend on some other utility (such as EZMAP) to have done it.

If CONPACKT calls SET, it always uses LNLG = 1, requesting a linear-linear mapping from the window to the viewport, and it positions the viewport and window as follows: The viewport is positioned as specified by the current values of the parameters 'VPL', 'VPR', 'VPB', 'VPT', and 'VPS'. The first four of these specify the position of a "viewport area", in which the viewport is to be centered and made as large as possible; the final one says how the shape of the viewport is to be determined. By default, the position of the window in the user coordinate system is determined by computing the minimum and maximum values of X and Y over all the points of the triangular mesh. The parameters 'WDL', 'WDR', 'WDB', and 'WDT' may be used to override this default behavior and specify the exact values to be used in the SET call to define the window.

If the triangular mesh represents data on the surface of a globe of radius 1, then to map the CONPACKT output onto an EZMAP background, one need only set 'MAP' to 1 and initialize EZMAP (which results in a call to the SPPS routine SET).

Special-Value Areas

CONPACK had a parameter 'SPV', whose value could be used in a data array in place of missing or unreliable data: any cell with this value at one or more of its vertices was essentially omitted from the rectangular grid. CONPACKT has no such parameter, but there are two ways to achieve the same effect:

There is no effective difference between triangles that are actually omitted from the mesh and triangles that are left in the mesh and marked as "blocked". The routine CTCLAM adds to an area map edges having omitted or blocked triangles on one side of them. The routines CTCLDM and CTCLDR may be made to draw edges having omitted or blocked triangles on one side of them (by giving a non-zero value to element "-1" of the parameter array 'CLU').

The Out-of-Range Parameter and Out-of-Range Areas

The parameter 'ORV', if non-zero, specifies an "out-of-range" value. This is only of use when the parameter 'MAP' is non-zero, specifying that coordinates are to be mapped by calling the subroutine CTMXYZ. The X coordinate returned by CTMXYZ may be set equal to 'ORV' to indicate that the mapped point is outside the range in which the mapping is defined.

A possible value for 'ORV', if it is to be set non-zero, is 1.E12, which has historically been returned by the EZMAP routines MAPTRN and MAPTRA to indicate a point which is outside the area depicted by a given map projection.

The union of all points for which CTMXYZ returns the out-of-range value constitutes a set of out-of-range areas. Contour lines are not traced in out-of-range areas (indeed, they cannot be). A binary-halving technique is used to extend contour lines to the very edge of such areas. The routine CTCLAM will attempt to generate and add to the area map a set of edges for such areas, and the routines CTCLDM and CTCLDR may be made to attempt to draw the edges of such areas (by giving a non-zero value to element "-2" of the parameter array 'CLU').

When contour lines are traced, if two consecutive points are out of range (in range), then the entire line segment connecting those two points is assumed to be out of range (in range). If the detail of the out-of-range areas is small enough, this assumption may cause errors. Giving the parameter 'PIC' a non-zero value will cause more points to be examined along each such line segment, thus curing the problem. For similar reasons, the algorithms used to trace the edge of the grid, the edges of special-value areas, and the edges of out-of-range areas may fail. Giving the parameter 'PIE' a non-zero value will cause these algorithms to use a finer grid, thus, again, curing the problem.

2D Smoothing of Contour Lines

Each of the routines which generates contour lines (CTCLAM, CTCLDM, and CTCLDR and the internal routine which positions labels along the lines using either the penalty scheme or the regular scheme) may be caused to smooth those lines, using cubic splines under tension. Giving the parameter 'T2D' a non-zero value causes this smoothing to be done. The absolute value of 'T2D' specifies the tension to be used; values near zero (.001, for example) yield approximately cubic splines (which may give very "loopy" curves), and large values (15., for example) yield nearly polygonal curves. Using very large values (50., for example) can cause overflow in the routines that do the smoothing.

Note that, since each contour line is smoothed separately, there is no way to absolutely ensure that it will not cause adjacent contour lines to cross each other; one must experiment with the tension to reduce the probability of that to a minimum. A reasonable value to start with is 2.5.

If 'T2D' is negative, smoothing is done prior to the coordinate mapping, if any, implied by the setting of the parameter 'MAP'. If 'T2D' is positive, smoothing is done after the mapping.

The parameter 'PIC' says how many points are to be interpolated between each pair of points defining the contour line, before smoothing. If 'PIC' is given a non-zero value when the 2D smoother is turned on, the effect is to constrain the smoothed curves to more closely approximate the original polygonal lines.

The parameter 'SSL' specifies the distance between points used to draw the smoothed contour lines. It is expressed as a fraction of the width of the window in the coordinate system in which the smoothing is being done.

Dash Package Use by CONPACKT

By default, CONPACKT uses routines from the old dash package called DASHCHAR. By using the appropriate flags on the "ncargf77" command line (or, if "ncargf77" is not being used, by specifying a different set of libraries on a link/load command), one can substitute one of the other members of the family (DASHSMTH or DASHSUPR).

One can also tell CONPACKT to use routines from a new dash package, called DASHPACK, which is intended to replace all of the older packages. One does this by giving the internal parameter 'DPU' a negative value in place of the default positive value. The new dash package offers many advantages over the older ones; for a full description of it, see the programmer document called "DASHPACK, A SOFTWARE PACKAGE FOR DRAWING DASHED LINES".

When DASHPACK is used, the following considerations apply:

Hachuring of Contour Lines

The parameters 'HCL', 'HCF', and 'HCS' may be used to request that contour lines be "hachured" in one of various ways. "Hachuring" is the process of drawing "hachures": short line segments perpendicular to the contour and usually drawn on the "downslope" side of it.

Note: When hachuring is activated, it may be necessary to increase the value of the internal parameter named 'RWC'. See the descriptions of the parameters 'HCF' and 'RWC' for more information.

Labels

Two or three different types of labels are written by a call to the routine CTLBDR: an informational label, high and low labels, and contour-line labels (the latter only if ABS('LLP') is greater than 1). Boxes surrounding these labels may be added to an area map by calling the routine CTLBAM; the purpose of this will probably be to prevent contour lines drawn by a subsequent call to the routine CTCLDM from passing through the labels or to prevent the label boxes from being filled or colored by a subsequent call to the routine ARSCAM (in the package AREAS).

When a constant field was detected by the initial call to CTMESH, a fourth type of label may be written by CONPACKT routines. In this case, a call to CTLBDR will write a constant-field label, warning of the situation, in place of the labels it would normally write. A call to CTLBAM will add the label box for the constant-field label to the area map, instead of the label boxes for the other labels. Calls to CTCLDR and CTCLDM which would normally draw contour lines will write the constant-field label instead.

The appearance of all of these labels may be determined in detail by setting parameters:

In all of the above parameter names, a suffixed 'A' means "angle", 'B' means "box flag", 'C' means "color index", 'L' means "line width", 'S' means "size of characters", 'T' means "text of label", 'W' means "white space width", 'X' means "X coordinate", and 'Y' means "Y coordinate". The 'O' in 'HLO' means "overlap", while the 'O' in 'LLO' means "orientation".

All labels are written by means of calls to the character-plotting routine PLCHHQ, in the package PLOTCHAR. The angle, in degrees, at which a label is written is determined by the value of the parameter 'xxA' (and, if it is a contour-line label, by the value of the parameter 'LLO'). The box flag 'xxB' determines whether or not, prior to writing the label, a box surrounding it is filled, and whether or not, after writing the label, the edge of the box is drawn. If the box is filled, it is done using the color index specified by the parameter 'LBC'; if the edge of the box is drawn, it is done using the color index, if any, chosen for the label itself, which is determined by the value of the parameter 'xxC'. The line width to be used in drawing the box is determined by the value of the parameter 'xxL'. The size (width) of the characters is determined by the value of the parameter 'xxS'. The text of the label is determined by the value of the parameter 'xxT'; usually, this string may contain embedded substrings of the form '$xxx$', which are to be replaced by the value of the quantity specified by the three-character mnemonic 'xxx'. The width of the "white space" to be left around the label (which defines the dimensions of the box around it) is determined by the value of the parameter 'xxW'.

Positioning of Labels on Contour Lines

CONPACKT provides three different ways to position labels on contour lines. The parameter 'LLP' is set by the user to choose one of the three.

When 'LLP' is 2 or 3, smoothing, if any, implied by a non-zero value of 'T2D' is suspended during placement of contour line labels. This saves a good deal of computer time and space and provides almost as good a set of labels. To leave smoothing turned on during placement of contour line labels, use 'LLP' = -2 or -3.

Details of the Penalty Scheme for Positioning Contour Line Labels

The penalty scheme is controlled by the ten parameters 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'PC6', 'PW1', 'PW2', 'PW3', and 'PW4'. (The rationale behind the names is that the first six parameters are characterized as "constants" and the remaining four as "weights" for the four terms in the penalty function.)

A point P on a contour line will be rejected as the center point of a label under any of the following conditions:

  1. If there are fewer than three points on the line.

  2. If the point P is too close to the center point of any other label on the current line, where the meaning of "too close" is defined by 'PC6'.

  3. If a label at the point P would extend outside the current viewport.

  4. If a label at the point P would overlap any previous label.

  5. If the estimated gradient of the field being contoured is too large at the point P, where "too large" is defined by the value of 'PC1', and 'PW1', the weight of the first term in the penalty function, is non-zero.

  6. If the estimated number of contour lines crossing a label at the point P is too large, where "too large" is defined by the value of 'PC2', and 'PW2', the weight of the second term in the penalty function, is non-zero.

  7. If the cumulative change in direction of the contour line within a circle covering a label at the point P is too large, where "too large" is defined by the value of 'PC3', and 'PW3', the weight of the third term in the penalty function, is non-zero.

The penalty function computed at each remaining point has the form

      PFUN = PW1 * GRAD / (GRAV+PC1*GRSD)            GRADIENT TERM
     +       PW2 * ENCB / PC2                        NUMBER-OF-CONTOURS TERM
     +       PW3 * CDIR / PC3                        CHANGE-IN-DIRECTION TERM
     +       PW4 * MIN (1-EXP(-((D(I)-PC4)/PC5)**2)) OPTIMUM-DISTANCE TERM
The first term of the penalty function becomes larger in high-gradient regions. GRAD is the estimated gradient at the point P, GRAV is the average gradient over the whole field being contoured, and GRSD is the standard deviation of the estimated gradients over the whole field. The parameter 'PC1' specifies how far from the norm gradients are allowed to wander, as a multiple of the standard deviation. Condition 5 above implies that, for points at which the penalty function is computed, either 'PW1' is zero or GRAD is less than or equal to GRAV+PC1*GRSD.

The second term of the penalty function becomes larger as ENCB, the estimated number of contour bands crossing a label at the point P, increases. The parameter 'PC2' specifies the largest number of crossing bands allowed. Condition 6 above implies that, for points at which the penalty function is computed, either 'PW2' is zero or ENCB is less than or equal to 'PC2'.

The third term of the penalty function becomes larger as CDIR, the cumulative change in direction of the contour line in a circular region centered at the point P and with a radius equal to half the larger dimension of the label, increases. The parameter 'PC3' specifies the largest such cumulative change allowed, in degrees. Condition 7 above implies that, for points at which the penalty function is computed, either 'PW3' is zero or CDIR is less than or equal to 'PC3'.

The fourth term of the penalty function becomes larger as the distance of the point P from the centers of all labels previously placed on other contour lines deviates from an optimum value specified by the user. D(I) represents the distance to the Ith such label center. The minimum is taken over all values of I. The parameter 'PC4' is the user-specified optimum distance, specified as a fraction of the width of the current viewport. If the point P is exactly 'PC4' units away from some previous label, then "MIN(1-EXP(...))" will have the value 0; otherwise, it will be non-zero. The parameter 'PC5' is the "folding distance", specified as a fraction of the width of the current viewport; as its value decreases, the function "1-EXP(...)" develops a sharper spike at D(I) = 'PC4'.

Thumbnail descriptions and default values of all the user-settable parameters are given below:

    'PC1'=1.   MULTIPLIER OF THE STANDARD DEVIATION OF THE GRADIENTS
    'PC2'=5.   MAXIMUM NUMBER OF CROSSING CONTOUR BANDS
    'PC3'=60.  MAXIMUM CUMULATIVE CHANGE IN DIRECTION OF THE CONTOUR LINE
    'PC4'=.05  OPTIMUM DISTANCE, AS A FRACTION OF THE WIDTH OF THE VIEWPORT
    'PC5'=.15  FOLDING DISTANCE, AS A FRACTION OF THE WIDTH OF THE VIEWPORT
    'PC6'=.30  MINIMUM DISTANCE BETWEEN LABELS ON THE SAME CONTOUR LINE, AS
                   A FRACTION OF THE WIDTH OF THE VIEWPORT
    'PW1'=2.   WEIGHT OF THE GRADIENT TERM
    'PW2'=0.   WEIGHT OF THE NUMBER-OF-CONTOURS TERM
    'PW3'=1.   WEIGHT OF THE CHANGE-IN-DIRECTION TERM
    'PW4'=1.   WEIGHT OF THE OPTIMUM-DISTANCE TERM

Choosing a Scale Factor

It is possible to specify a scale factor by which field values are to be divided before conversion to a character form for display as a numeric label.

The parameter 'SFS' says how the scale factor is to be chosen. If it is given a value greater than zero, that value is the desired scale factor. If 'SFS' is given a value less than or equal to zero, CONPACKT is directed to choose a scale factor to use, in one of five different ways. The value of the parameter 'SFU' may be retrieved by the user; it specifies the scale factor which has been selected for use.

The value of the scale factor may be displayed as a part of the informational label. This is done by embedding the substring '$SFU$' in the string which gives the value of the parameter 'ILT'.

The default value of 'SFS' is 1, which essentially specifies that no scale factor is to be used.

Constant Field Detection

The routine CTMESH checks for a field which is essentially constant. When such a field is found, the parameter 'CFF' is set non-zero; otherwise, 'CFF' is zeroed. The value of 'CFF' may be retrieved by the user.

When 'CFF' is non-zero, a call to one of the routines CTCLDM or CTCLDR will not cause any contour lines to be drawn; instead, the constant-field label will be written. The edge of the grid, the edges of special-value areas, and the edges of out-of-range areas may still be drawn.

Similarly, when 'CFF' is non-zero, a call to the routine CTLBDR will write the constant-field label instead of the labels which would normally be written, and a call to the routine CTLBAM will put the label box for the constant-field label into the area map instead of the label boxes for the normal set of labels.

Workspace Management

Many of the routines in CONPACKT require one or more workspace arrays, some of type REAL and some of type INTEGER. Some routines require more workspace under some conditions than under other conditions. (For example, when label positions are being chosen using either the regular scheme or the penalty scheme, space is required in both the real and integer workspaces and the amount required is larger for a complicated contour plot than it is for a simple one.) Building internal workspaces into a package causes problems; to enlarge such a workspace, the source code of the package must be modified, which is unacceptable. Providing each routine with arguments specifying all of the separate workspaces that it requires leads to complicated and error-prone code.

The workspace management scheme used in CONPACKT is as follows: The user defines one workspace array of type REAL and another of type INTEGER. In the call to CTMESH that initializes the drawing of a contour plot, these arrays appear as arguments (called RWRK and IWRK), together with arguments specifying their lengths (LRWK and LIWK). In subsequent calls to other CONPACKT routines which require workspaces, the same arrays appear as arguments, but the lengths do not. The CONPACKT routines cooperate in using these arrays in such a way as not to interfere with one another. Dynamic enlargement of one workspace at the expense of another becomes possible and the probability of running out of space is reduced.

In general, it is safest not to use the workspace arrays for other purposes between one call to a CONPACKT routine and the next (unless the next is to the routine CTMESH, which initializes the workspace pointers). At the moment, there is only one case in which the contents of the arrays are assumed to be preserved intact: If labels are being positioned using either the "regular" scheme or the "penalty" scheme, the list of label positions is created in the workspace arrays when it is first required and is assumed untouched thereafter. Other cases may arise as a result of further development of the package, however.

It is possible to find out how much space has been used in each of the workspace arrays. The parameters 'IWU' and 'RWU' are zeroed by a call to CTMESH and are updated thereafter to reflect maximum usage of space in the arrays. Thus, one might give the arrays large dimensions, create a typical contour plot, retrieve the values of 'IWU' and 'RWU' to see how much space was actually used, and then reduce the dimensions to more reasonable values.

Workspace usage by some routines cannot be predicted exactly. Upper bounds can be computed, but they may be rather large. For this reason, it may not be possible to absolutely ensure that enough workspace will be available for a given call. Therefore, there is a parameter, called 'WSO', which says what to do when a workspace overflow occurs. The four possibilities are as follows: to terminate after printing an error message, to continue running after printing an error message (this is the default), to just continue running without printing anything, or to do a recoverable-error call to SETER and then continue running without printing anything. Of course, in the latter two cases, incomplete plots may result. It is possible to find out whether or not a workspace overflow has occurred during a given call to a CONPACKT routine; this is done by retrieving the values of the parameters 'IWU' and 'RWU' and comparing them with the dimensions of the workspace arrays IWRK and RWRK.

The following information may be of value in attempting to estimate how much real and integer workspace will be required by each of the user-callable routines of CONPACKT. First, define the following quantities:

Then:

Routines not mentioned above use no workspace.

The Character-Width Multiplier

The parameter 'CWM' is used as a multiplier for all character widths and similar quantities. This makes it possible to scale all such quantities up and down simultaneously. The default value of 'CWM' is 1.

Searching for Highs and Lows

A point of the triangular mesh is defined to be the position of a high if and only if it is in the interior of the mesh and the data value there is greater than the data value at any other point of the mesh within a distance specified by the value of the internal parameter 'HLR'.

A similar definition is used for the position of a low.

Extended High/Low Search Algorithm

It is sometimes the case, in a data field being contoured (particularly if it has been packed and unpacked in such a way as to discretize the distribution of possible values in it), that the same value occurs at more than one mesh point in a small connected region A. If that value is greater than or equal to all the values in some extended region B around A, then it seems reasonable to place a high there, even though none of the mesh points in A satisfies the definition given above. Lows of a similar nature may also exist.

In January, 2002, CONPACK was modified to search for such highs and lows, but only if the user had set the value of the internal parameter 'HLE' non-zero. Such a feature was added to CONPACKT in July, 2004.

If the user has set the value of the internal parameter 'HLE' (which see) non-zero, an additional search will be performed, looking for connected regions in which the field value is constant. Each such region will then be tested to see if it ought to be considered a high or a low (or neither).

Setting 'HLE' equal to 1 will allow each small connected region A to be of any size. If 'HLE' is set larger than 1, it will set an upper limit on the size of each such region; for example, if 'HLE' = 4, then any such region containing more than 4 mesh points will be ignored.

The region B within which the field values must be smaller than or equal to (for a high) or greater than or equal to (for a low) the value at the mesh points of A will be just the union of the neighborhoods defined by the value of 'HLR' for all of the mesh points in A.

The algorithm to do this seems to be only moderately expensive to run, but its timing does depend on the nature of the user's contour field. The feature should probably not be used for fields having very large areas filled with a relative high or low. One problem arises that did not arise for CONPACK: if the high or low is spread out over a relatively large portion of the mesh where the mesh has a small radius of curvature, then the arithmetic mean of the X, Y, and Z coordinates of the points involved can yield a position for the high or low which is significantly off-mesh.

Color-Setting Philosophy

A number of the CONPACKT parameters specify the color of some object (like the informational label) or class of objects (like all contour lines for a given contour level). The default value of each such parameter is -1, which says that the color is to be determined by the current value of one of the GKS color indices. (The polyline color index is used for lines, the text color index for labels, and the fill area color index for filling label boxes.) If the value of the CONPACKT color-setting parameter for a given object is given a value greater than or equal to zero, it specifies the color index of the color in which the object is to be drawn. Before any object is drawn, the values of the GKS color indices affected are saved; after the object is drawn, the saved values are restored.

This structure allows the use of a tiered approach to color setting. If no color setting whatsoever is done, contour plots are drawn entirely in the colors specified by the applicable default values of the GKS color indices. If, on the other hand, prior to calling CONPACKT, one defines the color index "IC" (see the next section, "GKS Considerations") and then uses

      CALL GSPLCI (IC)
to change the GKS polyline color index, then all polylines drawn by CONPACKT change color. Similarly, one may use the statement

      CALL GSTXCI (IC)
to change the GKS text color index and the statement

    CALL GSFACI (IC)
to change the GKS fill area color index; the first will cause labels drawn by CONPACKT to change color and the second will cause label boxes filled by CONPACKT to change color.

If, in addition or instead, CONPACKT color-setting parameters are given values greater than or equal to zero, the objects or classes of objects to which those parameters apply are colored accordingly; these colors are used in preference to values preset by calls to GSPLCI, GSTXCI, or GSFACI.

A final opportunity to set color is provided by the user-supplied versions of the "change" routines, with names of the form CTCHxx; calls to GSPLCI, GSTXCI, and GSFACI may occur in such a routine and take precedence over color setting by any other means. Note that, if color is being set for drawing a label, then either or both of the polyline color index and the text color index may need to be set, depending on whether the labels are being drawn by calls to the GKS routine GPL (to draw polylines stroking out the characters) or by calls to the GKS routine GTX (to draw text). In particular, the routine PLCHHQ, in the package PLOTCHAR, which is called by CONPACKT to draw labels, may be directed by the user to draw high-quality characters, which are stroked, medium-quality characters, which are also stroked, or low-quality characters, which are drawn by GTX.

GKS Considerations

Certain assumptions are made by CONPACKT about the state of GKS, as follows:

(1) Like all the utilities in the NCAR graphics package, CONPACKT assumes that GKS has been opened and that the desired workstations have been opened and activated. The statement

      CALL OPNGKS
calls the SPPS routine OPNGKS, the GKS equivalent of which is

      CALL GOPKS (6,0)
      CALL GOPWK (1,2,1)
      CALL GACWK (1)
creating a single metacode workstation associated with FORTRAN unit 2.

Similarly, at the end of one's program, the workstations must be deactivated and closed and then GKS must be closed. The statement

    CALL CLSGKS
calls the SPPS routine CLSGKS, the GKS equivalent of which is

      CALL GDAWK (1)
      CALL GCLWK (1)
      CALL GCLKS
(2) It is assumed that the aspect source flags for various quantities are set to "individual". (The NCAR GKS package does this by default, but other packages may not.) To make sure that all the aspect source flags are set correctly, use the following code:

      DIMENSION IASF(13)
      ...
      DATA IASF / 13*1 /
      ...
      CALL GSASF (IASF)
(3) Color fill of label boxes is done by CONPACKT using calls to the GKS routine GFA; color fill of contour bands by the user will be done using similar calls. To get solid fill, rather than hollow fill, one must call a GKS routine to set the "fill area interior style":

      CALL GSFAIS (1)
(This is because the default "fill area interior style", as mandated by the GKS standard, is "hollow", rather than "solid".)

(4) Color-setting by CONPACKT is done by executing calls to the GKS routines GSPLCI, GSTXCI, and GSFACI, with user-defined color indices as arguments. The association of these color indices with colors on the workstations must have been defined previously by the user. This should be done by calling the GKS routine GSCR. The statement

      CALL GSCR (IW,IC,RC,GC,BC)
defines, for workstation IW, color index IC, with RGB components RC, GC, and BC. To be consistent with the SPPS routines OPNGKS and CLSGKS, use IW = 1. The value of IC may be any non-negative integer. By default, color index 0 is associated with the color black, which is defined by (RC,GC,BC) = (0.,0.,0.) and is used as the background color, while color index 1 is associated with the color white, which is defined by (RC,GC,BC) = (1.,1.,1.).

SUBROUTINES

All of the CONPACKT routines have six-character names beginning with the letters 'CT'. The user-callable ones are described in detail below.

CTBACK (RPNT,IEDG,ITRI,RWRK,IWRK)

This routine draws a background for a contour plot.

The initial version of CTBACK does very little. User feedback will be useful in determining what this routine is eventually made to do.

Usage

The routine CTBACK may be called at any time after the initialization call to CTMESH. All it does is draw the perimeter of the current viewport; it does this by calling PERIM, in the package GRIDAL.

Arguments

All five arguments are arrays used in a previous call to one of the routines CTMESH, CTMVRW, or CTMVIW.

RPNT (REAL array, dimensioned as specified in the last call to CTMESH, input) is the user's mesh-point array.

IEDG (INTEGER array, dimensioned as specified in the last call to CTMESH, input) is the user's edge array.

ITRI (INTEGER array, dimensioned as specified in the last call to CTMESH, input) is the user's triangle array.

RWRK (REAL array, dimensioned as specified in the last call to CTMESH or CTMVRW, input/output) is the current real workspace array.

IWRK (INTEGER array, dimensioned as specified in the last call to CTMESH or CTMVIW, input/output) is the current integer workspace array.

CTCHCF (IFLG)

This routine provides user control as a constant-field message is drawn.

Usage

The routine CTCHCF is not to be called by the user. It is called several times by CONPACKT while a constant-field label is being drawn. The default version of CTCHCF does nothing. A user-supplied replacement may change attributes such as color and line width (by calling the SPPS routine SETUSV or the appropriate GKS routines).

The text of the label being written may be retrieved by means of a "CALL CTGETC ('CTM',CVAL)". The text of the label may be changed by means of a "CALL CTSETC ('CTM',CVAL)"; this should only be done during a call with IFLG = 1 or 3 and, if it is done for one of those two values, it should also be done for the other.

When CTCHCF is called, the parameter 'DVA' will have been set to the value of the field; its value may be retrieved and used by CTCHCF.

Arguments

IFLG (INTEGER, input) is positive if an action is about to be taken, negative if an action has just been completed. The action in question is defined by the absolute value of IFLG, as follows:

CTCHCL (IFLG)

This routine provides user control as contour lines are drawn.

Usage

The routine CTCHCL is not to be called by the user. It is called by the CONPACKT routines CTCLDM and CTCLDR just before and just after the contour lines at each level are drawn. The default version of CTCHCL does nothing. A user-supplied replacement may change attributes such as color and line width (by calling the SPPS routine SETUSV or the appropriate GKS routines).

If the element of the parameter array 'CLU' corresponding to 'PAI' = -1 has been set non-zero to request the drawing of the edge of the grid, then CTCHCL will be called before and after that is done. Similarly, if the element of 'CLU' corresponding to 'PAI' = -2 has been set non-zero, then CTCHCL will be called before and after the drawing of the edges of the out-of-range areas.

When CTCHCL is called, the parameter 'PAI' will have been set to the index of the appropriate contour level (between 1 and 'NCL') or to one of the values -1 or -2. By retrieving the value of 'PAI', CTCHCL can find out what line is being drawn; also, a CTGETx call to retrieve an element of a parameter array like 'CLD' will automatically get the correct one for the line being drawn.

Arguments

IFLG (INTEGER, input) is +1 if a line is about to be drawn, -1 if a line has just been drawn.

CTCHHL (IFLG)

This routine provides user control as high and low labels are drawn.

Usage

The routine CTCHHL is not to be called by the user. It is called several times by CTLBAM and/or CTLBDR while each high or low label is positioned and drawn. The default version of CTCHHL does nothing. A user-supplied replacement may change attributes such as color and line width (by calling the SPPS routine SETUSV or the appropriate GKS routines).

The text of the label may be retrieved by means of a "CALL CTGETC ('CTM',CVAL)". The text of the label may be changed by means of a "CALL CTSETC ('CTM',CVAL)"; this should only be done during a call with IFLG = 1, 3, 5, or 7; if it is done for one of the two values 1 and 3, it should also be done for the other; if it is done for one of the two values 5 and 7, it should also be done for the other.

When CTCHHL is called, the parameter 'DVA' will have been set to the value of the high or low being labelled; its value may be retrieved and used by CTCHHL. Also, the internal parameters 'LBX' and 'LBY' will have been set to the X and Y coordinates of the center point of the label, in the current user coordinate system.

Arguments

IFLG (INTEGER, input) is positive if an action is about to be taken, negative if an action has just been completed. The action in question is defined by the absolute value of IFLG, as follows:

CTCHIL (IFLG)

This routine provides user control as the informational label is drawn.

Usage

The routine CTCHIL is not to be called by the user. It is called several times by CTLBAM and/or CTLBDR while the informational label is positioned and drawn. The default version of CTCHIL does nothing. A user-supplied replacement may change attributes such as color and line width (by calling the SPPS routine SETUSV or the appropriate GKS routines).

The text of the label may be retrieved by means of a "CALL CTGETC ('CTM',CVAL)". The text of the label may be changed by means of a "CALL CTSETC ('CTM',CVAL)"; this should only be done during a call with IFLG = 1 or IFLG = 3, and, if it is done for one of those values, it should be done for the other.

The internal parameters 'LBX' and 'LBY' will have been set to the X and Y coordinates of the center point of the label, in the current user coordinate system.

Arguments

IFLG (INTEGER, input) is positive if an action is about to be taken, negative if an action has just been completed. The action in question is defined by the absolute value of IFLG, as follows:

CTCHLL (IFLG)

This routine provides user control as contour line labels are drawn.

Usage

The routine CTCHLL is not to be called by the user. It is called several times by CTLBAM and/or CTLBDR while each contour-line label is positioned and drawn. The default version of CTCHLL does nothing. A user-supplied replacement may change attributes such as color and line width (by calling the SPPS routine SETUSV or the appropriate GKS routines).

The text of the label may be retrieved by means of a "CALL CTGETC ('CTM',CVAL)". The text of the label may be changed by means of a "CALL CTSETC ('CTM',CVAL)"; this should only be done during a call with IFLG = 1 or 3 and, if it is done for one of those values, it should be done for the other.

When CTCHLL is called, the parameter 'PAI' will have been set to the index of the appropriate contour level. Thus, parameters associated with that level may easily be retrieved by calls to CTGETx. Also, the parameter 'DVA' will have been set to the contour level value and the internal parameters 'LBX' and 'LBY' will have been set to the X and Y coordinates of the center point of the label, in the current user coordinate system.

Arguments

IFLG (INTEGER, input) is positive if an action is about to be taken, negative if an action has just been completed. The action in question is defined by the absolute value of IFLG, as follows:

CTCICA (RPNT,IEDG,ITRI,RWRK,IWRK, ... )

(The remaining arguments are ICRA,ICA1,ICAM, ICAN,XCPF, YCPF, XCQF, and YCQF.)

This routine creates a GKS cell array, using color indices determined by examining where the user's contours lie relative to the cell array. It ignores all triangles of the triangular mesh that are blocked, as specified by the current values of the blocking masks 'TBX' and 'TBA' and the blocking flag for each triangle. If the 3D package TDPACK is being used to project the triangular mesh into the plane, then the routine CTTDCA should probably be called instead, since it can be made to create a much better picture of the mesh than CTCICA can do.

Usage

The routine CTCICA may be called at any time after the initialization call to CTMESH. It is given arrays defining a triangular mesh of data, a real workspace array, an integer workspace array, an array in which a GKS cell array is to be created, and dimensioning and positioning information for the cell array. It creates the cell array by associating with each cell of it an area identifier (let's call it IAID) and then using IAID to generate a color index to store in that cell; there are three possible values of IAID:

Note that, to be sure of distinguishing all three possibilities from each other, it is necessary to give elements -1 and -2 of the parameter array 'AIA' values which are different from the value of any other element of 'AIA' or 'AIB'.

Note also that the "off-grid" value is used in cases in which it would not have been by the CONPACK routine CPCICA; in particular, it may be used for points inside mapped triangles of the triangular mesh having one or two vertices that map out-of-range.

Once an area identifier has been associated with a particular cell, that area identifier must be mapped into a color index to be stored as the value of the corresponding element of the cell array. The internal parameter 'CAF' determines how area identifiers are mapped into color indices, as follows:

The default value of 'CAF' is zero, so that CTSCAE is not called: if IAID is non-negative, it is used as the desired color index for a particular cell; otherwise, a zero is used.

Note that the routine CTCICA works differently from the routine CPCICA, in CONPACK, which allowed one to modify some elements of an existing cell array, while leaving others, in out-of-range areas or off-grid areas or special-value areas, unchanged. CPCICA could do that because the CONPACK routine CPMPXY implemented full inverse mappings, which the CONPACKT routine CTMXYZ does not do.

Arguments

The first five arguments are arrays used in a previous call to one of the routines CTMESH, CTMVRW, or CTMVIW.

RPNT (REAL array, dimensioned as specified in the last call to CTMESH, input) is the user's mesh-point array.

IEDG (INTEGER array, dimensioned as specified in the last call to CTMESH, input) is the user's edge array.

ITRI (INTEGER array, dimensioned as specified in the last call to CTMESH, input) is the user's triangle array.

RWRK (REAL array, dimensioned as specified in the last call to CTMESH or CTMVRW, input/output) is the current real workspace array.

IWRK (INTEGER array, dimensioned as specified in the last call to CTMESH or CTMVIW, input/output) is the current integer workspace array.

ICRA (INTEGER array, dimensioned ICA1 by "n", where "n" is greater than or equal to the value of the argument ICAN, input/output) is the user's cell array.

ICA1 (INTEGER, input) is the first dimension of the FORTRAN array ICRA, which contains the user's cell array.

ICAM (INTEGER, input) is the first dimension of the user's cell array.

ICAN (INTEGER, input) is the second dimension of the user's cell array.

XCPF and YCPF (REAL, input) are the coordinates, in the fractional coordinate system, of a point P, which is the point at that corner of the rectangular area into which the cell array maps that is also a corner of the cell with indices (1,1).

XCQF and YCQF (REAL, input) are the coordinates, in the fractional coordinate system, of a point Q, which is the point at that corner of the rectangular area into which the cell array maps that is also a corner of the cell with indices (ICAM,ICAN).

CTCLAM (RPNT,IEDG,ITRI,RWRK,IWRK,IAMA)

This routine adds contour lines to an area map as part of the process of drawing a solid-fill contour plot. It ignores triangles of the triangular mesh that are blocked, as specified by the current values of the blocking masks 'TBX' and 'TBA' and the blocking flag for each triangle.

Usage

The routine CTCLAM may be called at any time after the initialization call to CTMESH to add contour lines generated from the data defined on the triangular mesh to the area map in the array IAMA. The area map must previously have been initialized by a call to the routine ARINAM, in the package AREAS.

The contour lines added to the area map are as specified by the first 'NCL' elements of the parameter arrays 'CLV', 'AIA', and 'AIB'. If 'NCL' is zero, CTPKCL is called to generate these values. The contour levels defined by the first 'NCL' elements of the parameter array 'CLV' are then examined. If a given contour level is associated with a non-zero value of 'AIA' and/or a non-zero value of 'AIB', the contour lines at that contour level are added to the area map. If there is an associated non-zero value of 'AIA', it is used as the area identifier for the area "above" the line (where field values are greater than they are along the line); otherwise, a zero is used. If there is an associated non-zero value of 'AIB', it is used as the area identifier for the area "below" the line (where field values are less than they are along the line); otherwise, a zero is used. Note that a given contour level may occur more than once in the internal parameter array 'CLV', but there must be at most one non-zero value of 'AIA' and at most one non-zero value of 'AIB' associated with it; otherwise, an error exit occurs. If the parameter 'T2D' has a non-zero value, the contour lines are smoothed, using cubic splines under tension.

Other types of lines are also added to the area map by CTCLAM: the edge of the current viewport and possibly a set of vertical lines within the viewport; the edge of the grid; and the edges of out-of-range areas, if any. The area identifier for the outside of the viewport is always -1. Elements of the parameter array 'AIA' for 'PAI' = -1 and -2 may be used to specify the area identifiers to be used for the outside of the grid and the inside of an out-of-range area, respectively; the default values of the first is 0 and the default value of the second is -1. Area identifiers for all other sides of these edges are determined from the area-identifier information given for the contour levels.

The lines added to the area map are put into two edge groups, one with group identifier 'GIC' and another with group identifier 'GIS'. (The edge of the viewport may actually be added twice, once to each group.) The edge group 'GIC' is the more important of the two; it really defines the division of the plane into contour bands, out-of-range areas, and outside-the-grid areas. The edge group 'GIS' only receives the edge of the viewport and a collection of vertical lines; its object is to break up the areas defined by the other edge group into vertical pieces, creating simpler polygons on devices that might not handle more complicated ones. Whether or not edge group 'GIS' is created is under control of the user; for more information, see the descriptions of the parameters 'GIS' and 'NVS'.

Lines are added to the area map in the following order: the edge of the viewport and the vertical lines within it, edges of the out-of-range areas, if any; the edge of the grid; and the contour lines, in order of increasing contour level.

If, during the last call to CTMESH, the data being contoured were found to be essentially constant, then no contour lines are added to the area map; the other lines are added, however.

Arguments

The first five arguments are arrays used in a previous call to the one of the routines CTMESH, CTMVRW, or CTMVIW.

RPNT (REAL array, dimensioned as specified in the last call to CTMESH, input) is the user's mesh-point array.

IEDG (INTEGER array, dimensioned as specified in the last call to CTMESH, input) is the user's edge array.

ITRI (INTEGER array, dimensioned as specified in the last call to CTMESH, input) is the user's triangle array.

RWRK (REAL array, dimensioned as specified in the last call to CTMESH or CTMVRW, input/output) is the current real workspace array.

IWRK (INTEGER array, dimensioned as specified in the last call to CTMESH or CTMVIW, input/output) is the current integer workspace array.

IAMA (INTEGER array, dimensioned as specified in a call to ARINAM, in the package AREAS) is the array containing the area map to which contour lines are to be added.

CTCLDM (RPNT,IEDG,ITRI,RWRK,IWRK,IAMA,RTPL)

This routine draws contour lines masked by an existing area map. The object of this may be simply to avoid drawing contour lines through label boxes, but the routine may be used for more complicated tasks, like limiting the drawing of contour lines to the ocean areas on an EZMAP background. It ignores triangles of the triangular mesh that are blocked, as specified by the current values of the blocking masks 'TBX' and 'TBA' and the blocking flag for each triangle.

Usage

The routine CTCLDM may be called at any time after the initialization call to CTMESH to draw contour lines masked by an existing area map. Actually, CTCLDM does not draw the lines; it generates them, masks them against a user-specified area map, and generates a set of calls to a user-specified routine (one call for each polyline resulting from the masking process). Each such polyline lies entirely within precisely one of the areas defined by the area map. The user routine may use the information provided by its arguments, describing the area the polyline is in, to decide whether or not to draw the polyline.

The contour lines generated are those specified by the first 'NCL' elements of the parameter arrays 'CLV' and 'CLU'. If 'NCL' is zero, CTPKCL is called to generate these values. Each element of 'CLV' specifies a contour level and the corresponding element of 'CLU' specifies whether or not contour lines are to be generated at that level. If the parameter 'T2D' has a non-zero value, the contour lines are smoothed, using cubic splines under tension.

If the eleme