PROGRAM CTCBAY
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
PARAMETER (IERF=6,LUNI=2,IWTY=1 ,IWID=1) ! NCGM
C PARAMETER (IERF=6,LUNI=2,IWTY=8 ,IWID=1) ! X Windows
C PARAMETER (IERF=6,LUNI=2,IWTY=20,IWID=1) ! PostScript
C PARAMETER (IERF=6,LUNI=2,IWTY=11,IWID=1) ! PDF, Portrait
C PARAMETER (IERF=6,LUNI=2,IWTY=12,IWID=1) ! PDF, Landscape
C
C The object of this program is to produce a set of plots illustrating
C the use of a triangular mesh derived from a grid on the surface of
C the globe provided by Tom Gross.
C
C Up to four frames are drawn at each selected viewpoint:
C
C 1) A frame showing the triangular mesh on the globe. This frame
C is drawn only if the parameter IMSH is non-zero.
C
C 2) A frame showing simple contours on the globe. This frame is
C drawn only if the parameter ICON is non-zero.
C
C 3) A frame showing color-filled contour bands on the globe, drawn
C using filled areas. This frame is drawn only if the parameter
C ICOL is non-zero.
C
C 4) A frame showing color-filled contour bands on the globe, drawn
C using a cell array. This frame is drawn only if the parameter
C ICAP is non-zero.
C
C Define the parameter that says whether or not the triangular mesh is
C to be drawn (frame 1):
C
C PARAMETER (IMSH=0) ! triangular mesh not drawn
PARAMETER (IMSH=1) ! triangular mesh drawn
C
C Define the parameter that says whether or not to draw simple contours
C (frame 2):
C
C PARAMETER (ICON=0) ! contours not drawn
PARAMETER (ICON=1) ! contours drawn
C
C Define the parameter that says whether or not to draw color-filled
C contours (frame 3):
C
C PARAMETER (ICOL=0) ! no color fill done
PARAMETER (ICOL=1) ! color fill done
C
C Define the parameter that says whether or not to draw a cell array
C plot (frame 4):
C
C PARAMETER (ICAP=0) ! cell array plot not drawn
PARAMETER (ICAP=1) ! cell array plot drawn
C
C To represent the triangular mesh, we use three singly-dimensioned
C arrays: RPNT holds points, IEDG holds edges, and ITRI holds triangles.
C The elements of each array form "nodes" having lengths as follows:
C
PARAMETER (LOPN=5) ! length of a point node
C
C The five elements of a point node are
C
C 1. the X coordinate of the point;
C 2. the Y coordinate of the point;
C 3. the Z coordinate of the point;
C 4. the field value at the point;
C 5. any additional value desired by the user.
C
PARAMETER (LOEN=5) ! length of an edge node
C
C The five elements of an edge node are
C
C 1. the base index, in RPNT, of point 1 of the edge;
C 2. the base index, in RPNT, of point 2 of the edge;
C 3. the index, in ITRI, of the pointer to the edge in the triangle to
C the left of the edge (-1 if there is no triangle to the left);
C 4. the index, in ITRI, of the pointer to the edge in the triangle to
C the right of the edge (-1 if there is no triangle to the right);
C 5. a utility flag for use by algorithms that scan the structure.
C
C The "left" and "right" sides of an edge are defined as they would be
C by an observer standing on the globe at point 1 of the edge, looking
C toward point 2 of the edge. It is possible, if there are "holes" in
C the mesh, that there will be no triangle to the left or to the right
C of an edge, but there must be a triangle on one side or the other.
C
PARAMETER (LOTN=4) ! length of a triangle node
C
C The four elements of a triangle node are
C
C 1. the base index, in IEDG, of edge 1 of the triangle;
C 2. the base index, in IEDG, of edge 2 of the triangle;
C 3. the base index, in IEDG, of edge 3 of the triangle;
C 4. a flag set non-zero to block use of the triangle, effectively
C removing it from the mesh. Use the ISSCP grid and play with
C the setting of the parameter ISCP (which see, below) to get
C examples of the use of this feature.
C
C The "base index" of a point node, an edge node, or a triangle node is
C always a multiple of the length of the node, to which can be added an
C offset to get the index of a particular element of the node. For
C example, if I is the base index of a triangle of interest, ITRI(I+1)
C is its first element (the base index of its first edge). Similarly,
C IEDG(ITRI(I+1)+2) is the base index of the second point of the first
C edge of the triangle with base index I, and RPNT(IEDG(ITRI(I+1)+2)+3)
C is the third (Z) coordinate of the second point of the first edge of
C the triangle with base index I.
C
C It is the pointers from the edge nodes back to the triangle nodes that
C allow CONPACKT to navigate the mesh, moving from triangle to triangle
C as it follows a contour line, but these pointers are tricky to define:
C if IPTE is the base index of an edge node and IEDG(IPTE+3) is zero or
C more, saying that there is a triangle to the left of the edge, then
C IEDG(IPTE+3) is the actual index of that element of the triangle node
C that points to the edge node; i.e., ITRI(IDGE(IPTE+3))=IPTE. The base
C index of the triangle node defining that triangle is IPTT, where
C IPTT=LOTN*((IEDG(IPTE+3)-1)/LOTN), and the index of the pointer to
C the edge within the triangle node is IPTI=IEDG(IPTE+3)-IPTT, so that
C IEDG(IPTT+IPTI)=IPTE. Similar comments apply to element 4 of an edge
C node, which points into the triangle node defining the triangle to the
C right of the edge.
C
C The numbers of points, edges, and triangles in the triangular mesh are
C approximately as follows (obtained from a run in which the numbers
C were set larger than this):
C
PARAMETER (MNOP=10000)
PARAMETER (MNOE=23000)
PARAMETER (MNOT=16000)
C
C Once we know how many points, edges, and triangles we're going to use
C (at most), we can set parameters defining the space reserved for the
C triangular mesh:
C
PARAMETER (MPNT=MNOP*LOPN) ! space for points
PARAMETER (MEDG=MNOE*LOEN) ! space for edges
PARAMETER (MTRI=MNOT*LOTN) ! space for triangles
C
C Declare the arrays to hold the point nodes, edge nodes, and triangle
C nodes defining the triangular mesh.
C
DIMENSION RPNT(MPNT),IEDG(MEDG),ITRI(MTRI)
C
C Declare sort arrays to be used in GTCBAY to keep track of where points
C and edges were put in the structure defining the triangular mesh.
C
DIMENSION IPPP(2,MNOP),IPPE(2,MNOE)
C
C Declare real and integer workspaces needed by CONPACKT.
C
PARAMETER (LRWK=10000,LIWK=1000)
C
DIMENSION RWRK(LRWK),IWRK(LIWK)
C
C Declare the area map array needed to do solid fill.
C
PARAMETER (LAMA=400000)
C
DIMENSION IAMA(LAMA)
C
C Declare workspace arrays to be used in calls to ARSCAM.
C
PARAMETER (NCRA=LAMA/10,NGPS=2)
C
DIMENSION XCRA(NCRA),YCRA(NCRA),IAAI(NGPS),IAGI(NGPS)
C
C Declare arrays in which to generate a cell array picture of the
C data on the triangular mesh.
C
C PARAMETER (ICAM=1024,ICAN=1024)
PARAMETER (ICAM= 512,ICAN= 512)
C
DIMENSION ICRA(ICAM,ICAN)
C
C Declare external the routine that draws masked contour lines.
C
EXTERNAL DRWMCL
C
C Declare external the routine that does color fill of contour bands.
C
EXTERNAL DCFOCB
C
C Define a common block in which to keep track of the maximum space used
C in the arrays XCRA and YCRA.
C
COMMON /COMONA/ MAXN
C
C Define the out-of-range flag.
C
DATA OORV / 1.E12 /
C
C Define the tension on the splines to be used in smoothing contours.
C
DATA T2DS / 0.0 / ! smoothing off
C DATA T2DS / 2.5 / ! smoothing on
C
C Define the distance between points on smoothed contour lines.
C
DATA RSSL / .002 /
C
C Define the amount of real workspace to be used in drawing contours.
C
DATA IRWC / 500 /
C
C Define the label-positioning flag.
C
C DATA ILLP / 0 / ! no labels
C DATA ILLP / 1 / ! dash-package writes labels
DATA ILLP / 2 / ! regular scheme
C DATA ILLP / 3 / ! penalty scheme
C
C Define the high/low search radius.
C
DATA HLSR / .075 /
C
C Define the high/low label overlap flag.
C
DATA IHLO / 11 /
C
C Define the hachuring flag, hachure length, and hachure spacing.
C
DATA IHCF,HCHL,HCHS / 0 , +.004 , .010 / ! off
C DATA IHCF,HCHL,HCHS / +1 , -.004 , .020 / ! on, all, uphill
C
C Define a constant to convert from radians to degrees.
C
DATA RTOD / 57.2957795130823 /
C
C Read data and generate the required triangular mesh.
C
PRINT * , ' '
PRINT * , 'CREATING TRIANGULAR MESH:'
C
CALL GTCBAY (RPNT,MPNT,NPNT,LOPN, ! point list
+ IEDG,MEDG,NEDG,LOEN, ! edge list
+ ITRI,MTRI,NTRI,LOTN, ! triangle list
+ IPPP,IPPE,CLAT,CLON)
C
C Print the number of points, edges, and triangles.
C
PRINT * , ' NUMBER OF POINTS: ',NPNT/LOPN
PRINT * , ' NUMBER OF EDGES: ',NEDG/LOEN
PRINT * , ' NUMBER OF TRIANGLES: ',NTRI/LOTN
C
C Write the contents of the point list, the edge list, and the triangle
C list to "fort.11" in a readable form.
C
C WRITE (11,'(''P'',I8,5F10.4)')
C + (I,RPNT(I+1),RPNT(I+2),RPNT(I+3),RPNT(I+4),RPNT(I+5),
C + I=0,NPNT-LOPN,LOPN)
C WRITE (11,'(''E'',I8,5I10)')
C + (I,IEDG(I+1),IEDG(I+2),IEDG(I+3),IEDG(I+4),IEDG(I+5),
C + I=0,NEDG-LOEN,LOEN)
C WRITE (11,'(''T'',I8,4I10)')
C + (I,ITRI(I+1),ITRI(I+2),ITRI(I+3),ITRI(I+4),
C + I=0,NTRI-LOTN,LOTN)
C
C Open GKS.
C
PRINT * , 'OPENING AND INITIALIZING GKS'
C
CALL GOPKS (IERF,0)
CALL GOPWK (IWID,LUNI,IWTY)
CALL GACWK (IWID)
C
C Turn off the clipping indicator.
C
CALL GSCLIP (0)
C
C Define a basic set of colors.
C
CALL GSCR (IWID, 0,1.,1.,1.) ! white (background)
CALL GSCR (IWID, 1,0.,0.,0.) ! black (foreground)
CALL GSCR (IWID, 2,1.,1.,0.) ! yellow
CALL GSCR (IWID, 3,1.,0.,1.) ! magenta
CALL GSCR (IWID, 4,1.,0.,0.) ! red
CALL GSCR (IWID, 5,0.,1.,1.) ! cyan
CALL GSCR (IWID, 6,0.,1.,0.) ! green
CALL GSCR (IWID, 7,0.,0.,1.) ! blue
CALL GSCR (IWID, 8,.5,.5,.5) ! darker light gray
CALL GSCR (IWID, 9,.8,.8,.8) ! lighter light gray
CALL GSCR (IWID,10,.3,.3,0.) ! dark yellow
CALL GSCR (IWID,11,.3,.3,.3) ! dark gray
CALL GSCR (IWID,12,.5,.5,.5) ! medium gray
CALL GSCR (IWID,13,.8,.5,.5) ! light red - geographic lines
CALL GSCR (IWID,14,.5,.5,.8) ! light blue - lat/lon lines
C
C Define 100 colors, associated with color indices 151 through 250, to
C be used for color-filled contour bands and in cell arrays, ranging
C from blue to red.
C
CALL DFCLRS (IWID,151,250,0.,0.,1.,1.,0.,0.)
C
C Set parameters in the utilities.
C
PRINT * , 'SETTING PARAMETERS IN CONPACKT, EZMAP, AND PLOTCHAR'
C
C Set the mapping flag.
C
CALL CTSETI ('MAP - MAPPING FLAG',1)
C
C Set the out-of-range flag value.
C
CALL CTSETR ('ORV - OUT-OF-RANGE VALUE',OORV)
C
C Turn on the drawing of the mesh edge and set the area identifier for
C areas outside the mesh.
C
CALL CTSETI ('PAI - PARAMETER ARRAY INDEX',-1)
CALL CTSETI ('CLU - CONTOUR LEVEL USE FLAG',1)
C CALL CTSETI ('AIA - AREA IDENTIFIER FOR AREA',1001)
C
C Set the area identifier for areas in "out-of-range" areas.
C
C CALL CTSETI ('PAI',-2)
C CALL CTSETI ('AIA - AREA IDENTIFIER FOR AREA',1002)
C
C Set the 2D smoother flag.
C
CALL CTSETR ('T2D - TENSION ON 2D SPLINES',T2DS)
C
C Set the distance between points on smoothed lines.
C
CALL CTSETR ('SSL - SMOOTHED SEGMENT LENGTH',RSSL)
C
C Set the amount of real workspace to be used in drawing contours.
C
CALL CTSETI ('RWC - REAL WORKSPACE FOR CONTOURS',IRWC)
C
C Set the label-positioning flag.
C
CALL CTSETI ('LLP - LINE LABEL POSITIONING FLAG',ILLP)
C
C Set the high/low search radius.
C
CALL CTSETR ('HLR - HIGH/LOW SEARCH RADIUS',HLSR)
C
C Set the high/low label overlap flag.
C
CALL CTSETI ('HLO - HIGH/LOW LABEL OVERLAP FLAG',IHLO)
C
C Set the hachuring flag, hachure length, and hachure spacing.
C
CALL CTSETI ('HCF - HACHURING FLAG',IHCF)
CALL CTSETR ('HCL - HACHURE LENGTH',HCHL)
CALL CTSETR ('HCS - HACHURE SPACING',HCHS)
C
C Set the cell array flag.
C
CALL CTSETI ('CAF - CELL ARRAY FLAG',-1)
C
C Tell CONPACKT not to do its own call to SET, since EZMAP will have
C done it.
C
CALL CTSETI ('SET - DO-SET-CALL FLAG', 0)
C
C Move the informational label up a little.
C
CALL CTSETR ('ILY - INFORMATIONAL LABEL Y POSITION',-.005)
C
C Tell EZMAP not to draw the perimeter.
C
CALL MPSETI ('PE',0)
C
C Put lat/lon lines at 1-degree intervals.
C
CALL MPSETI ('GR',1)
C
C Tell EZMAP to use solid lat/lon lines.
C
C CALL MPSETI ('DA',65535)
C
C Tell PLOTCHAR to use one of the filled fonts and to outline each
C character.
C
CALL PCSETI ('FN',25)
CALL PCSETI ('OF',1)
C
C Tell PLOTCHAR to expect a character other than a colon as the
C function-control signal character.
C
CALL PCSETC ('FC','|')
C
C Loop through the desired viewing angles.
C
DO 104 IDIR=1,1 ! one (OR) views of the globe
C
PRINT * , ' '
PRINT * , 'VIEW FROM DIRECTION NUMBER: ',IDIR
C
C Tell EZMAP what projection to use and what its limits are. (The ".2"
C added to CLAT is a heuristic modification to compensate for the fact
C that the center of the mesh is pulled down a bit due to the way in
C which I compute it.)
C
IF (IDIR.EQ.1) THEN
CALL MAPROJ ('OR',CLAT+.2,CLON,0.)
END IF
C
CALL MAPSET ('AN',1.6,1.6,1.6,1.6)
C
C Initialize EZMAP.
C
CALL MAPINT
C
C If the triangular mesh is to be drawn, do it.
C
IF (IMSH.NE.0) THEN
C
PRINT * , 'DRAWING TRIANGULAR MESH'
C
C Triangular mesh (with edges between blocked triangles somewhat fainter
C than the rest):
C
DO 101 IPTE=0,NEDG-LOEN,LOEN
C
IFLL=0
C
IF (IEDG(IPTE+3).GE.0) THEN
IF (ITRI(LOTN*((IEDG(IPTE+3)-1)/LOTN)+4).EQ.0) IFLL=1
END IF
C
IFLR=0
C
IF (IEDG(IPTE+4).GE.0) THEN
IF (ITRI(LOTN*((IEDG(IPTE+4)-1)/LOTN)+4).EQ.0) IFLR=1
END IF
C
IF (IFLL.NE.0.OR.IFLR.NE.0) THEN
CALL PLOTIT (0,0,2)
CALL GSPLCI (8)
ELSE
CALL PLOTIT (0,0,2)
CALL GSPLCI (9)
END IF
C
ALAT=RTOD*ASIN(RPNT(IEDG(IPTE+1)+3))
C
IF (RPNT(IEDG(IPTE+1)+1).EQ.0..AND.
+ RPNT(IEDG(IPTE+1)+2).EQ.0.) THEN
ALON=0.
ELSE
ALON=RTOD*ATAN2(RPNT(IEDG(IPTE+1)+2),
+ RPNT(IEDG(IPTE+1)+1))
END IF
C
BLAT=RTOD*ASIN(RPNT(IEDG(IPTE+2)+3))
C
IF (RPNT(IEDG(IPTE+2)+1).EQ.0..AND.
+ RPNT(IEDG(IPTE+2)+2).EQ.0.) THEN
BLON=0.
ELSE
BLON=RTOD*ATAN2(RPNT(IEDG(IPTE+2)+2),
+ RPNT(IEDG(IPTE+2)+1))
END IF
C
CALL DRSGCR (ALAT,ALON,BLAT,BLON)
C
101 CONTINUE
C
CALL PLOTIT (0,0,2)
CALL GSPLCI (13)
CALL MAPGRD
CALL PLOTIT (0,0,2)
CALL GSPLCI (14)
CALL MDLNDR ('Earth..1',1)
C
C Label the first frame.
C
CALL PLCHHQ (CFUX(.03),CFUY(.946),'CHESAPEAKE BAY',
+ .024,0.,-1.)
CALL PLCHHQ (CFUX(.03),CFUY(.904),'FINITE ELEMENT MESH',
+ .018,0.,-1.)
C
CALL PLCHHQ (CFUX(.97),CFUY(.950),'TRIANGULAR MESH',
+ .012,0.,1.)
C
C Advance the frame.
C
CALL FRAME
C
END IF
C
C If a frame showing simple contours is to be drawn, do it next (adding
C an overlay of lat/lon lines and continental outlines in light gray).
C
IF (ICON.NE.0) THEN
C
PRINT * , 'DRAWING PLOT SHOWING SIMPLE CONTOURS'
C
C Initialize CONPACKT.
C
PRINT * , 'CALLING CTMESH'
C
CALL CTMESH (RPNT,NPNT,LOPN, ! point list
+ IEDG,NEDG,LOEN, ! edge list
+ ITRI,NTRI,LOTN, ! triangle list
+ RWRK,LRWK, ! real workspace
+ IWRK,LIWK) ! integer workspace
C
CALL PLOTIT (0,0,2)
CALL GSPLCI (1)
C
C Proceed as implied by the setting of the label-positioning flag.
C
IF (ABS(ILLP).EQ.1) THEN
C
C Draw the contour lines with labels generated by the dash package.
C
PRINT * , 'CALLING CTCLDR'
CALL CTCLDR (RPNT,IEDG,ITRI,RWRK,IWRK)
C
C Add the informational and high/low labels.
C
PRINT * , 'CALLING CTLBDR'
CALL CTLBDR (RPNT,IEDG,ITRI,RWRK,IWRK)
C
ELSE IF (ABS(ILLP).GT.1) THEN
C
C Create an area map for masking of labels.
C
MAXN=0
C
PRINT * , 'CALLING ARINAM'
CALL ARINAM (IAMA,LAMA)
C
PRINT * , 'CALLING CTLBAM'
CALL CTLBAM (RPNT,IEDG,ITRI,RWRK,IWRK,IAMA)
C
C Draw the contour lines masked by the area map.
C
PRINT * , 'CALLING CTCLDM'
CALL CTCLDM (RPNT,IEDG,ITRI,RWRK,IWRK,IAMA,DRWMCL)
C
PRINT * , 'AREA MAP SPACE REQUIRED: ',
+ IAMA(1)-IAMA(6)+IAMA(5)+1
C
PRINT * , 'NUMBER OF POINTS IN LONGEST LINE:',MAXN
C
C Draw all the labels.
C
PRINT * , 'CALLING CTLBDR'
CALL CTLBDR (RPNT,IEDG,ITRI,RWRK,IWRK)
C
END IF
C
CALL CTGETI ('IWU - INTEGER WORKSPACE USED',IIWU)
PRINT * , 'INTEGER WORKSPACE REQUIRED: ',IIWU
C
CALL CTGETI ('RWU - REAL WORKSPACE USED',IRWU)
PRINT * , 'REAL WORKSPACE REQUIRED: ',IRWU
C
C Add lat/lon lines and continental outlines.
C
CALL PLOTIT (0,0,2)
CALL GSPLCI (13)
CALL MAPGRD
CALL PLOTIT (0,0,2)
CALL GSPLCI (14)
CALL MDLNDR ('Earth..1',1)
C
C Label the second frame.
C
CALL PLCHHQ (CFUX(.03),CFUY(.946),'CHESAPEAKE BAY',
+ .024,0.,-1.)
CALL PLCHHQ (CFUX(.03),CFUY(.904),'FINITE ELEMENT MESH',
+ .018,0.,-1.)
C
CALL PLCHHQ (CFUX(.97),CFUY(.950),'SIMPLE CONTOURS ON',
+ .012,0.,1.)
CALL PLCHHQ (CFUX(.97),CFUY(.926),'TRIANGULAR MESH',
+ .012,0.,1.)
C
CALL PLCHHQ (CFUX(.97),CFUY(.904),'Data: temperatures at',
+ .010,0.,1.)
CALL PLCHHQ (CFUX(.97),CFUY(.884),'time 1, level 1.',
+ .010,0.,1.)
C
C Advance the frame.
C
CALL FRAME
C
END IF
C
C If a frame showing color-filled contours is to be drawn, do it next.
C
IF (ICOL.NE.0) THEN
C
PRINT * , 'DRAWING PLOT SHOWING COLOR-FILLED CONTOURS'
C
PRINT * , 'CALLING CTMESH'
C
CALL CTMESH (RPNT,NPNT,LOPN, ! point list
+ IEDG,NEDG,LOEN, ! edge list
+ ITRI,NTRI,LOTN, ! triangle list
+ RWRK,LRWK, ! real workspace
+ IWRK,LIWK) ! integer workspace
C
MAXN=0
C
PRINT * , 'CALLING CTPKCL'
CALL CTPKCL (RPNT,IEDG,ITRI,RWRK,IWRK)
C
PRINT * , 'CALLING ARINAM'
CALL ARINAM (IAMA,LAMA)
C
PRINT * , 'CALLING CTCLAM'
CALL CTCLAM (RPNT,IEDG,ITRI,RWRK,IWRK,IAMA)
C
C PRINT * , 'CALLING CTLBAM'
C CALL CTLBAM (RPNT,IEDG,ITRI,RWRK,IWRK,IAMA)
C
PRINT * , 'CALLING ARSCAM'
CALL ARSCAM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,DCFOCB)
C
PRINT * , 'SPACE REQUIRED IN AREA MAP: ',
+ IAMA(1)-IAMA(6)+IAMA(5)+1
C
PRINT * , 'NUMBER OF POINTS IN LARGEST AREA:',MAXN
C
CALL CTGETI ('IWU - INTEGER WORKSPACE USED',IIWU)
PRINT * , 'INTEGER WORKSPACE REQUIRED: ',IIWU
C
CALL CTGETI ('RWU - REAL WORKSPACE USED',IRWU)
PRINT * , 'REAL WORKSPACE REQUIRED: ',IRWU
C
CALL GSFACI (1)
CALL PLOTIT (0,0,2)
CALL GSPLCI (13)
CALL MAPGRD
CALL PLOTIT (0,0,2)
CALL GSPLCI (14)
CALL MDLNDR ('Earth..1',1)
C
C Label the third frame.
C
CALL PLCHHQ (CFUX(.03),CFUY(.946),'CHESAPEAKE BAY',
+ .024,0.,-1.)
CALL PLCHHQ (CFUX(.03),CFUY(.904),'FINITE ELEMENT MESH',
+ .018,0.,-1.)
C
CALL PLCHHQ (CFUX(.97),CFUY(.950),'COLORED CONTOUR BANDS ON'
+,.012,0.,1.)
CALL PLCHHQ (CFUX(.97),CFUY(.926),'TRIANGULAR MESH',
+ .012,0.,1.)
C
CALL PLCHHQ (CFUX(.97),CFUY(.904),'Data: temperatures at',
+ .010,0.,1.)
CALL PLCHHQ (CFUX(.97),CFUY(.884),'time 1, level 1.',
+ .010,0.,1.)
C
C Advance the frame.
C
CALL FRAME
C
END IF
C
C If the flag for it is set, do a cell array plot.
C
IF (ICAP.NE.0) THEN
C
PRINT * , 'DRAWING CELL-ARRAY PLOT OF DATA VALUES'
C
PRINT * , 'CALLING CTMESH'
C
CALL CTMESH (RPNT,NPNT,LOPN, ! point list
+ IEDG,NEDG,LOEN, ! edge list
+ ITRI,NTRI,LOTN, ! triangle list
+ RWRK,LRWK, ! real workspace
+ IWRK,LIWK) ! integer workspace
C
CALL GETSET (XVPL,XVPR,YVPB,YVPT,XWDL,XWDR,YWDB,YWDT,LNLG)
C
PRINT * , 'CALLING CTCICA'
CALL CTCICA (RPNT,IEDG,ITRI,RWRK,IWRK,ICRA,ICAM,ICAM,ICAN,
+ XVPL,YVPB,XVPR,YVPT)
C
PRINT * , 'CALLING GCA'
CALL GCA (XWDL,YWDB,XWDR,YWDT,ICAM,ICAN,1,1,ICAM,ICAN,ICRA)
C
CALL CTGETI ('IWU - INTEGER WORKSPACE USED',IIWU)
PRINT * , 'INTEGER WORKSPACE REQUIRED: ',IIWU
C
CALL CTGETI ('RWU - REAL WORKSPACE USED',IRWU)
PRINT * , 'REAL WORKSPACE REQUIRED: ',IRWU
C
CALL PLOTIT (0,0,2)
CALL GSPLCI (13)
CALL MAPGRD
CALL PLOTIT (0,0,2)
CALL GSPLCI (14)
CALL MDLNDR ('Earth..1',1)
C
C Label the fourth frame.
C
CALL PLCHHQ (CFUX(.03),CFUY(.946),'CHESAPEAKE BAY',
+ .024,0.,-1.)
CALL PLCHHQ (CFUX(.03),CFUY(.904),'FINITE ELEMENT MESH',
+ .018,0.,-1.)
C
CALL PLCHHQ (CFUX(.97),CFUY(.950),'CELL ARRAY DERIVED FROM',
+.012,0.,1.)
CALL PLCHHQ (CFUX(.97),CFUY(.926),'TRIANGULAR MESH',
+ .012,0.,1.)
C
CALL PLCHHQ (CFUX(.97),CFUY(.904),'Data: temperatures at',
+ .010,0.,1.)
CALL PLCHHQ (CFUX(.97),CFUY(.884),'time 1, level 1.',
+ .010,0.,1.)
C
C Advance the frame.
C
CALL FRAME
C
END IF
C
104 CONTINUE
C
C Close GKS.
C
CALL GDAWK (IWID)
CALL GCLWK (IWID)
CALL GCLKS
C
C Done.
C
STOP
C
END
SUBROUTINE DFCLRS (IWID,IOFC,IOLC,REDF,GRNF,BLUF,REDL,GRNL,BLUL)
C
C This routine defines color indices IOFC through IOLC on workstation
C IWID by interpolating values from REDF/GRNF/BLUF to REDL/GRNL/BLUL.
C
DO 101 I=IOFC,IOLC
P=REAL(IOLC-I)/REAL(IOLC-IOFC)
Q=REAL(I-IOFC)/REAL(IOLC-IOFC)
CALL GSCR (IWID,I,P*REDF+Q*REDL,P*GRNF+Q*GRNL,P*BLUF+Q*BLUL)
101 CONTINUE
C
C Done.
C
RETURN
C
END
SUBROUTINE DRWMCL (XCRA,YCRA,NCRA,IAAI,IAGI,NGPS)
C
DIMENSION XCRA(*),YCRA(*),IAAI(*),IAGI(*)
C
C This routine draws the curve defined by the points (XCRA(I),YCRA(I)),
C for I = 1 to NCRA, if and only if none of the area identifiers for the
C area containing the polyline are negative. It calls either CURVE or
C CURVED to do the drawing, depending on the value of the internal
C parameter 'DPU'.
C
C It keeps track of the maximum value used for NCRA in a common block.
C
COMMON /COMONA/ MAXN
C
MAXN=MAX(MAXN,NCRA)
C
C Retrieve the value of the internal parameter 'DPU'.
C
CALL CTGETI ('DPU - DASH PACKAGE USED',IDUF)
C
C Turn on drawing.
C
IDRW=1
C
C If any area identifier is negative, turn off drawing.
C
DO 101 I=1,NGPS
IF (IAAI(I).LT.0) IDRW=0
101 CONTINUE
C
C If drawing is turned on, draw the polyline.
C
IF (IDRW.NE.0) THEN
IF (IDUF.EQ.0) THEN
CALL CURVE (XCRA,YCRA,NCRA)
ELSE IF (IDUF.LT.0) THEN
CALL DPCURV (XCRA,YCRA,NCRA)
ELSE
CALL CURVED (XCRA,YCRA,NCRA)
END IF
END IF
C
C Done.
C
RETURN
C
END
SUBROUTINE DCFOCB (XCRA,YCRA,NCRA,IAAI,IAGI,NGPS)
C
C This routine fills the area defined by the points (XCRA(I),YCRA(I)),
C for I = 1 to NCRA, if and only if none of the area identifiers for
C the area are negative. The color used is determined from the area
C identifier of the area relative to group 3; it is assumed that 100
C colors are defined having color indices 151 through 250.
C
DIMENSION XCRA(*),YCRA(*),IAAI(*),IAGI(*)
C
C It keeps track of the maximum value used for NCRA in a common block.
C
COMMON /COMONA/ MAXN
C
MAXN=MAX(MAXN,NCRA)
C
C Retrieve the number of contour levels being used.
C
CALL CTGETI ('NCL - NUMBER OF CONTOUR LEVELS',NOCL)
C
C If the number of contour levels is non-zero and the area has more
C than two points, fill it.
C
IF (NOCL.NE.0.AND.NCRA.GT.2) THEN
C
IAI3=-1
C
DO 101 I=1,NGPS
IF (IAGI(I).EQ.3) IAI3=IAAI(I)
101 CONTINUE
C
IF (IAI3.GE.1.AND.IAI3.LE.NOCL+1) THEN
CALL GSFACI (151+INT(((REAL(IAI3)-.5)/REAL(NOCL+1))*100.))
CALL GFA (NCRA,XCRA,YCRA)
ELSE IF (IAI3.EQ.1001) THEN
CALL GSFACI (2)
CALL GFA (NCRA,XCRA,YCRA)
ELSE IF (IAI3.EQ.1002) THEN
CALL GSFACI (3)
CALL GFA (NCRA,XCRA,YCRA)
END IF
C
END IF
C
C Done.
C
RETURN
C
END
SUBROUTINE GTCBAY (RPNT,MPNT,NPNT,LOPN,
+ IEDG,MEDG,NEDG,LOEN,
+ ITRI,MTRI,NTRI,LOTN,
+ IPPP,IPPE,CLAT,CLON)
C
DIMENSION RPNT(MPNT),IEDG(MEDG),ITRI(MTRI),IPPP(2,*),IPPE(2,*)
C
C Construct a triangular mesh using data from Tom Gross, using a routine
C (CTTMTL) that makes it easy to create an arbitrary triangular mesh.
C The "quicksorts" used by this method are very inefficient for objects
C that are partially ordered, so, as the triangles of the mesh are
C generated, they are stored in a triangle buffer from which they can
C be dumped in random order by calls to the routine CTTMTL.
C
C Declare the array to use as the buffer for the triangles, so that the
C order in which they are processed can be somewhat randomized. Up to
C a point, making this buffer larger will result in more randomization
C of the triangles and speed up the process. MBUF is the number of
C triangles that will fit in the buffer at once, and KBUF is the number
C of those that are dumped whenever the buffer is found to be full; it
C turns out that it's better to dump more than 1 at a time (I suppose
C because the call and loop set-up time for CTTMTL are non-trivial, so
C it's better to dump a few triangles while you're there). So ... use
C MBUF > KBUF > 1.
C
PARAMETER (MBUF=5021,KBUF=173,DTOR=.017453292519943)
C
DIMENSION TBUF(12,MBUF)
C
C Include definitions that "NetCDF" needs (now commented out, because
C the data are being read from an ASCII file, instead).
C
C include 'netcdf.inc'
C
C Declare arrays to receive lat/lon coordinate data, contour field data,
C and triangle vertex indices from Tom's NetCDF data file.
C
C DIMENSION TLAT(7258),TLON(7258),TDAT(7258,11,25),IELE(3,13044)
C
C Declare arrays to receive lat/lon coordinate data, contour field data,
C and triangle vertex indices from the ASCII file made from Tom's data.
C
DIMENSION TLAT(7258),TLON(7258),TDAT(7258 ),IELE(3,13044)
C
C The data for this example were originally read from a user's "NetCDF"
C file. The code used for the purpose follows, but has been commented
C out, as the data are now read from an ASCII file. This gets around
C certain procedural problems in running the example from "ncargex".
C
C Open the "NetCDF" file.
C
C ISTA=NF_OPEN('ctcbay.nc',0,NCID)
C
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_OPEN: ',ISTA
C STOP
C END IF
C
C Read the array of latitudes of point nodes.
C
C ISTA=NF_INQ_VARID(NCID,'lat',IVID)
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_INQ_VARID: ',ISTA
C STOP
C END IF
C ISTA=NF_GET_VAR_REAL(NCID,IVID,TLAT)
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_GET_VAR_REAL: ',ISTA
C STOP
C END IF
C
C Read the array of longitudes of point nodes.
C
C ISTA=NF_INQ_VARID(NCID,'lon',IVID)
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_INQ_VARID: ',ISTA
C STOP
C END IF
C ISTA=NF_GET_VAR_REAL(NCID,IVID,TLON)
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_GET_VAR_REAL: ',ISTA
C STOP
C END IF
C
C Read an array of data at the point nodes.
C
C ISTA=NF_INQ_VARID(NCID,'temp',IVID) ! temperatures
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_INQ_VARID: ',ISTA
C STOP
C END IF
C ISTA=NF_GET_VAR_REAL(NCID,IVID,TDAT)
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_GET_VAR_REAL: ',ISTA
C STOP
C END IF
C
C Read the array of pointers to point triplets.
C
C ISTA=NF_INQ_VARID(NCID,'ele',IVID)
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_INQ_VARID: ',ISTA
C STOP
C END IF
C ISTA=NF_GET_VAR_INT(NCID,IVID,IELE)
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_GET_VAR_INT: ',ISTA
C STOP
C END IF
C
C Close the "NetCDF" file.
C
C ISTA=NF_CLOSE(NCID)
C IF (ISTA.NE.NF_NOERR) THEN
C PRINT * , 'ERROR RETURN FROM NF_CLOSE: ',ISTA
C STOP
C END IF
C
C Write the data to unit 11. (These statements are commented out; they
C were used to create a "fort.11" which was then renamed "ctcbay.dat".)
C
C WRITE (11,'(5E16.8)') TLAT
C WRITE (11,'(5E16.8)') TLON
C WRITE (11,'(5E16.8)') (TDAT(I,1,1),I=1,7258)
C WRITE (11, '(20I4)') IELE
C
C Read the required data from the ASCII file "ctcbay.dat".
C
OPEN (11,FILE='ctcbay.dat',STATUS='OLD',FORM='FORMATTED')
C
READ (11,'(5E16.8)') TLAT
READ (11,'(5E16.8)') TLON
READ (11,'(5E16.8)') TDAT
READ (11, '(20I4)') IELE
C
CLOSE (11)
C
C Build structures forming the triangular mesh. First, initialize the
C variables used to keep track of items in the sort arrays for points
C and edges, in the triangle buffer, and in the triangle list.
C
MPPP=MPNT/LOPN ! point sorting
NPPP=0
C
MPPE=MEDG/LOEN ! edge sorting
NPPE=0
C
NBUF=0 ! triangle buffer
C
NTRI=0 ! triangle list
C
C Loop through the triangles derived from Tom's grid.
C
DO 101 I=1,13044
C
C If the triangle buffer is almost full, dump a few randomly-chosen
C triangles from it to make room for two new ones.
C
IF (NBUF.GE.MBUF) THEN
CALL CTTMTL (KBUF,TBUF,MBUF,NBUF,
+ IPPP,MPPP,NPPP,
+ IPPE,MPPE,NPPE,
+ RPNT,MPNT,NPNT,LOPN,
+ IEDG,MEDG,NEDG,LOEN,
+ ITRI,MTRI,NTRI,LOTN)
END IF
C
C Put a new triangle in the triangle buffer.
C
NBUF=NBUF+1
C
TBUF( 1,NBUF)=COS(DTOR*TLAT(IELE(1,I)))*
+ COS(DTOR*TLON(IELE(1,I)))
TBUF( 2,NBUF)=COS(DTOR*TLAT(IELE(1,I)))*
+ SIN(DTOR*TLON(IELE(1,I)))
TBUF( 3,NBUF)=SIN(DTOR*TLAT(IELE(1,I)))
C TBUF( 4,NBUF)=TDAT(IELE(1,I),1,1) ! (data from NetCDF file)
TBUF( 4,NBUF)=TDAT(IELE(1,I) ) ! (data from ASCII file)
C
TBUF( 5,NBUF)=COS(DTOR*TLAT(IELE(2,I)))*
+ COS(DTOR*TLON(IELE(2,I)))
TBUF( 6,NBUF)=COS(DTOR*TLAT(IELE(2,I)))*
+ SIN(DTOR*TLON(IELE(2,I)))
TBUF( 7,NBUF)=SIN(DTOR*TLAT(IELE(2,I)))
C TBUF( 8,NBUF)=TDAT(IELE(2,I),1,1) ! (data from NetCDF file)
TBUF( 8,NBUF)=TDAT(IELE(2,I) ) ! (data from ASCII file)
C
TBUF( 9,NBUF)=COS(DTOR*TLAT(IELE(3,I)))*
+ COS(DTOR*TLON(IELE(3,I)))
TBUF(10,NBUF)=COS(DTOR*TLAT(IELE(3,I)))*
+ SIN(DTOR*TLON(IELE(3,I)))
TBUF(11,NBUF)=SIN(DTOR*TLAT(IELE(3,I)))
C TBUF(12,NBUF)=TDAT(IELE(3,I),1,1) ! (data from NetCDF file)
TBUF(12,NBUF)=TDAT(IELE(3,I) ) ! (data from ASCII file)
C
101 CONTINUE
C
C Output all triangles remaining in the buffer.
C
IF (NBUF.NE.0) THEN
CALL CTTMTL (NBUF,TBUF,MBUF,NBUF,
+ IPPP,MPPP,NPPP,
+ IPPE,MPPE,NPPE,
+ RPNT,MPNT,NPNT,LOPN,
+ IEDG,MEDG,NEDG,LOEN,
+ ITRI,MTRI,NTRI,LOTN)
END IF
C
C Set values telling the caller how many points and edges were created.
C
NPNT=NPPP*LOPN
NEDG=NPPE*LOEN
C
C Compute average X, Y, and Z values, from which we can get the
C approximate center point of the mesh.
C
XSUM=0.
YSUM=0.
ZSUM=0.
C
DO 102 I=0,NPNT-LOPN,LOPN
XSUM=XSUM+RPNT(I+1)
YSUM=YSUM+RPNT(I+2)
ZSUM=ZSUM+RPNT(I+3)
102 CONTINUE
C
XAVG=XSUM/REAL(NPPP)
YAVG=YSUM/REAL(NPPP)
ZAVG=ZSUM/REAL(NPPP)
C
DNOM=SQRT(XAVG*XAVG+YAVG*YAVG+ZAVG*ZAVG)
C
XAVG=XAVG/DNOM
YAVG=YAVG/DNOM
ZAVG=ZAVG/DNOM
C
C Return the latitude and longitude of the approximate center point of
C the mesh on the globe.
C
CLAT=ASIN(ZAVG)/DTOR
CLON=ATAN2(YAVG,XAVG)/DTOR
C
C Done.
C
RETURN
C
END
SUBROUTINE DRSGCR (ALAT,ALON,BLAT,BLON)
C
C (DRSGCR = DRaw Shortest Great Circle Route)
C
C This routine draws the shortest great circle route joining two points
C on the globe.
C
PARAMETER (MPTS=181,SIZE=1.) ! MPTS = INT(180./SIZE) + 1
C
DIMENSION QLAT(MPTS),QLON(MPTS)
C
NPTS=MAX(1,MIN(MPTS,INT(ADSGCR(ALAT,ALON,BLAT,BLON)/SIZE)))
C
CALL MAPGCI (ALAT,ALON,BLAT,BLON,NPTS,QLAT,QLON)
C
CALL MAPIT (ALAT,ALON,0)
C
DO 101 I=1,NPTS
CALL MAPIT (QLAT(I),QLON(I),1)
101 CONTINUE
C
CALL MAPIT (BLAT,BLON,1)
C
CALL MAPIQ
C
RETURN
C
END
FUNCTION ADSGCR (ALAT,ALON,BLAT,BLON)
C
C (ADSGCR = Angle in Degrees along Shortest Great Circle Route)
C
C This function returns the shortest great circle distance, in degrees,
C between two points, A and B, on the surface of the globe.
C
DATA DTOR / .017453292519943 /
DATA RTOD / 57.2957795130823 /
C
XCOA=COS(DTOR*ALAT)*COS(DTOR*ALON)
YCOA=COS(DTOR*ALAT)*SIN(DTOR*ALON)
ZCOA=SIN(DTOR*ALAT)
C
XCOB=COS(DTOR*BLAT)*COS(DTOR*BLON)
YCOB=COS(DTOR*BLAT)*SIN(DTOR*BLON)
ZCOB=SIN(DTOR*BLAT)
C
ADSGCR=2.*RTOD*ASIN(SQRT((XCOA-XCOB)**2+
+ (YCOA-YCOB)**2+
+ (ZCOA-ZCOB)**2)/2.)
C
RETURN
C
END
SUBROUTINE CTSCAE (ICRA,ICA1,ICAM,ICAN,XCPF,YCPF,XCQF,YCQF,
+ IND1,IND2,ICAF,IAID)
DIMENSION ICRA(ICA1,*)
C
C This routine is called by CTCICA when the internal parameter 'CAF' is
C given a negative value. Each call is intended to create a particular
C element in the user's cell array. The arguments are as follows:
C
C ICRA is the user's cell array.
C
C ICA1 is the first dimension of the FORTRAN array ICRA.
C
C ICAM and ICAN are the first and second dimensions of the cell array
C stored in ICRA.
C
C (XCPF,YCPF) is the point at that corner of the rectangular area
C into which the cell array maps that corresponds to the cell (1,1).
C The coordinates are given in the fractional coordinate system (unlike
C what is required in a call to GCA, in which the coordinates of the
C point P are in the world coordinate system).
C
C (XCQF,YCQF) is the point at that corner of the rectangular area into
C which the cell array maps that corresponds to the cell (ICAM,ICAN).
C The coordinates are given in the fractional coordinate system (unlike
C what is required in a call to GCA, in which the coordinates of the
C point Q are in the world coordinate system).
C
C IND1 is the 1st index of the cell that is to be updated.
C
C IND2 is the 2nd index of the cell that is to be updated.
C
C ICAF is the current value of the internal parameter 'CAF'. This
C value will always be an integer which is less than zero (because
C when 'CAF' is zero or greater, this routine is not called).
C
C IAID is the area identifier associated with the cell. It will have
C been given one of the values from the internal parameter array 'AIA'
C (the one for 'PAI' = -2 if the cell lies in an out-of-range area, the
C one for 'PAI' = -1 if the cell lies off the data grid, or the one for
C some value of 'PAI' between 1 and 'NCL' if the cell lies on the data
C grid). The value zero may occur if the cell falls in an out-of-range
C area and the value of 'AIA' for 'PAI' = -2 is 0 or if the cell lies
C off the data grid and the value of 'AIA' for 'PAI' = -1 is 0, or if
C the cell falls on the data grid, but no contour level below the cell
C has a non-zero 'AIA' and no contour level above the cell has a
C non-zero 'AIB'. Note that, if the values of 'AIA' for 'PAI' = -1
C and -2 are given non-zero values, IAID can only be given a zero
C value in one way.
C
C The default behavior of CTSCAE is as follows: If the area identifier
C is non-negative, it is treated as a color index, to be stored in the
C appropriate cell in the cell array; but if the area identifier is
C negative, a zero is stored for the color index. The user may supply
C a version of CTSCAE that does something different; it may simply map
C the area identifiers into color indices or it may somehow modify the
C existing cell array element to incorporate the information provided
C by the area identifier.
C
C ICRA(IND1,IND2)=MAX(0,IAID)
C
C What follows is not the default behavior; instead, it is the behavior
C expected by many of the example programs for CONPACKT.
C
CALL CTGETI ('NCL - NUMBER OF CONTOUR LEVELS',NOCL)
C
IF (IAID.GE.1.AND.IAID.LE.NOCL+1) THEN
ICRA(IND1,IND2)=151+INT(((REAL(IAID)-.5)/REAL(NOCL+1))*100.)
ELSE IF (IAID.EQ.1001) THEN
ICRA(IND1,IND2)=2
ELSE IF (IAID.EQ.1002) THEN
ICRA(IND1,IND2)=3
END IF
C
RETURN
C
END