

      PROGRAM CTORCA
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 to represent an ORCA grid on the surface
C of the globe.  The triangular mesh is created from the rectangular
C grid by calling the CONPACKT routine CTTMRG.
C
C Selected frames are drawn for each of four different viewpoints.  If
C (CLAT,CLON) is the approximate position of the "center point" of the
C mesh (which is defined somewhat arbitrarily for certain meshes), the
C four viewpoints used are as follows: (CLAT+45,CLON), (CLAT-45,CLON),
C (CLAT+45,CLON+180), and (CLAT-45,CLON+180).
C
C At each of the four viewpoints, up to six different frames are drawn:
C
C   1) A frame that attempts to show how the original rectangular
C      grid is mapped onto the globe.  A little reference image of
C      the rectangular grid appears in the lower left corner of the
C      frame and a large image of the grid wrapped around the globe
C      appears in the middle of the frame.  Each image is colored using
C      a 10x10 rectangular array of different colors and, on each, the
C      corner points, the midpoints of the edges and the midpoint of
C      the whole grid are marked.  This frame is drawn only if the
C      parameter IMSH = 1, 3, 5, or 7.
C
C   2) A frame showing the original rectangular grid on the globe.
C      This frame is drawn only if the parameter IMSH = 2, 3, 6, or 7.
C
C   3) A frame showing the triangular mesh on the globe.  This frame
C      is drawn only if the parameter IMSH = 4, 5, 6, or 7.
C
C   4) A frame showing simple contours on the globe.  This frame is
C      drawn only if the parameter ICON is non-zero.
C
C   5) 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   6) 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 which frames showing the rectangular
C grid and/or the triangular mesh are to be drawn (frames 1 through 3):
C
C       PARAMETER (IMSH=0)  !  none of the three
C       PARAMETER (IMSH=1)  !  colored rectangular grid only
C       PARAMETER (IMSH=2)  !  simple rectangular grid only
C       PARAMETER (IMSH=3)  !  colored and simple rectangular grids
C       PARAMETER (IMSH=4)  !  triangular only
C       PARAMETER (IMSH=5)  !  colored rectangular and triangular
C       PARAMETER (IMSH=6)  !  simple rectangular and triangular
        PARAMETER (IMSH=7)  !  all three
C
C Define the parameter that says whether or not to draw simple contours
C (frame 4):
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 5):
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 6):
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 maximum numbers of points, edges, and triangles in a triangular
C mesh can be computed in the following manner:  Since the triangular
C mesh is created from a rectangular one by splitting each rectangular
C cell in half, we can declare the largest values we expect to use for
C the dimensions of that mesh:
C
        PARAMETER (IDIM=181,JDIM=148,IDM1=IDIM-1,JDM1=JDIM-1)
C
C and then compute from those values the maximum number of points,
C edges, and triangles that the triangular mesh arrays will need to
C hold.  The computed values will be exactly as required if no points
C or edges of the rectangular grid are repeated in it; if there are
C repeating points or edges, space for slightly fewer points and edges
C will be needed:
C
        PARAMETER (MNOP=IDIM*JDIM)
        PARAMETER (MNOE=3*IDM1*JDM1+IDM1+JDM1)
        PARAMETER (MNOT=2*IDM1*JDM1)
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 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=200000)
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 arrays to hold color components for colors to be used at the
C four corners of an illustrative drawing of the rectangular grid.
C
        DIMENSION CCLL(3),CCLR(3),CCUL(3),CCUR(3)
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 the colors to be used in the lower left, lower right, upper
C left, and upper right corners of the rectangular grid when a drawing
C of it is made for the purpose of illustrating how it is wrapped around
C the globe.
C
        DATA CCLL / 0.2 , 0.2 , 0.2 /
        DATA CCLR / 1.0 , 0.0 , 1.0 /
        DATA CCUL / 0.0 , 1.0 , 1.0 /
        DATA CCUR / 1.0 , 1.0 , 0.0 /
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 GTORCA (RPNT,MPNT,NPNT,LOPN,  !  point list
     +               IEDG,MEDG,NEDG,LOEN,  !  edge list
     +               ITRI,MTRI,NTRI,LOTN,  !  triangle list
     +               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 51 through 150, to
C be used in making an illustrative drawing of the rectangular grid.
C
        CALL DCFGRD (IWID,51,10,10,CCLL,CCLR,CCUL,CCUR)
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)
        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 Tell EZMAP to use solid lat/lon lines.
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 four different viewing angles.
C
        DO 104 IDIR=1,4  !  four different (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.
C
          IF      (IDIR.EQ.1) THEN
            CALL MAPROJ ('OR',CLAT+45.,CLON,      0.)
          ELSE IF (IDIR.EQ.2) THEN
            CALL MAPROJ ('OR',CLAT-45.,CLON,      0.)
          ELSE IF (IDIR.EQ.3) THEN
            CALL MAPROJ ('OR',CLAT+45.,CLON+180., 0.)
          ELSE IF (IDIR.EQ.4) THEN
            CALL MAPROJ ('OR',CLAT-45.,CLON+180., 0.)
          END IF
C
          CALL MAPSET ('MA',0.,0.,0.,0.)
C
C Initialize EZMAP.
C
          CALL MAPINT
C
C If the rectangular grid and/or the triangular mesh are to be drawn,
C do it.
C
          IF (IMSH.NE.0) THEN
C
            PRINT * , 'DRAWING MESHES'
C
C Colored rectangular grid:
C
            IF (MOD(IMSH,2).NE.0) THEN
C
              PRINT * , 'DRAWING COLORED RECTANGULAR GRID'
C
              CALL DCORCA
C
C Label the first frame.
C
              CALL PLCHHQ (CFUX(.03),CFUY(.946),'ORCA GRID',.024,0.,-1.)
C
              IF      (IDIR.EQ.1) THEN
                CALL PLCHHQ (CFUX(.47),CFUY(.82),
     +                       'The top edge of the grid maps to a',
     +                                                       .010,0.,1.)
                CALL PLCHHQ (CFUX(.47),CFUY(.80),
     +                       'line passing through the North Pole.',
     +                                                       .010,0.,1.)
                CALL PLCHHQ (CFUX(.64),CFUY(.75),
     +                       'Portions of the grid over Asia and',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.64),CFUY(.73),
     +                       'Africa are shifted to provide better',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.64),CFUY(.71),
     +                       'coverage of various bodies of water.',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.64),CFUY(.69),
     +                       '(See viewpoint 3 for a better view.)',
     +                                                      .010,0.,-1.)
              ELSE IF (IDIR.EQ.2) THEN
                CALL PLCHHQ (CFUX(.53),CFUY(.30),
     +                       'The bottom edge of the grid maps',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.53),CFUY(.28),
     +                       'to a circle around the South Pole.',
     +                                                      .010,0.,-1.)
              ELSE IF (IDIR.EQ.3) THEN
                CALL PLCHHQ (CFUX(.53),CFUY(.77),
     +                       'The top edge of the grid maps to a',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.53),CFUY(.75),
     +                       'line passing through the North Pole.',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.40),CFUY(.45),
     +                       'Portions of the grid over Asia and',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.40),CFUY(.43),
     +                       'Africa are shifted to provide better',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.40),CFUY(.41),
     +                       'coverage of various bodies of water.',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.53),CFUY(.27),
     +                       'Left and right edges of',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.53),CFUY(.25),
     +                       'the grid map to a "seam"',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.53),CFUY(.23),
     +                       'on one of the meridians.',
     +                                                      .010,0.,-1.)
              ELSE IF (IDIR.EQ.4) THEN
                CALL PLCHHQ (CFUX(.40),CFUY(.82),
     +                       'Portions of the grid over Asia and',
     +                                                       .010,0.,1.)
                CALL PLCHHQ (CFUX(.40),CFUY(.80),
     +                       'Africa are shifted to provide better',
     +                                                       .010,0.,1.)
                CALL PLCHHQ (CFUX(.40),CFUY(.78),
     +                       'coverage of various bodies of water.',
     +                                                       .010,0.,1.)
                CALL PLCHHQ (CFUX(.40),CFUY(.76),
     +                       '(See viewpoint 3 for a better view.)',
     +                                                       .010,0.,1.)
                CALL PLCHHQ (CFUX(.47),CFUY(.37),
     +                       'The bottom edge of the grid maps',
     +                                                       .010,0.,1.)
                CALL PLCHHQ (CFUX(.47),CFUY(.35),
     +                       'to a circle around the South Pole.',
     +                                                       .010,0.,1.)
                CALL PLCHHQ (CFUX(.53),CFUY(.62),
     +                       'Left and right edges of',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.53),CFUY(.60),
     +                       'the grid map to a "seam"',
     +                                                      .010,0.,-1.)
                CALL PLCHHQ (CFUX(.53),CFUY(.58),
     +                       'on one of the meridians.',
     +                                                      .010,0.,-1.)
              END IF
C
              IF      (IDIR.EQ.1) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 1',
     +                                                      .012,0.,-1.)
              ELSE IF (IDIR.EQ.2) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 2',
     +                                                      .012,0.,-1.)
              ELSE IF (IDIR.EQ.3) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 3',
     +                                                      .012,0.,-1.)
              ELSE IF (IDIR.EQ.4) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 4',
     +                                                      .012,0.,-1.)
              END IF
C
              CALL PLCHHQ (CFUX(.97),CFUY(.950),'MAPPING OF THE GRID',
     +                                                       .012,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.926),'ONTO THE GLOBE',
     +                                                       .012,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.904),'(Color swatches, bold'
     +                                                      ,.010,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.884),'edges, and big dots',
     +                                                       .010,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.864),'are all intended',
     +                                                       .010,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.844),'to clarify the',
     +                                                       .010,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.824),'mapping.)',
     +                                                       .010,0.,1.)
C
              CALL PLCHHQ (CFUX(.97),CFUY(.060),'Shorelines are blue.',
     +                                                       .012,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.036),'Parallels/meridians are
     + red.',.012,0.,1.)
C
C Advance the frame.
C
              CALL FRAME
C
            END IF
C
C Simple rectangular grid:
C
            IF (MOD(IMSH/2,2).NE.0) THEN
C
              PRINT * , 'DRAWING SIMPLE RECTANGULAR GRID'
C
              CALL DSORCA
C
C Label the second frame.
C
              CALL PLCHHQ (CFUX(.03),CFUY(.946),'ORCA GRID',.024,0.,-1.)
C
              IF (     IDIR.EQ.1) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 1',
     +                                                      .012,0.,-1.)
              ELSE IF (IDIR.EQ.2) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 2',
     +                                                      .012,0.,-1.)
              ELSE IF (IDIR.EQ.3) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 3',
     +                                                      .012,0.,-1.)
              ELSE IF (IDIR.EQ.4) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 4',
     +                                                      .012,0.,-1.)
              END IF
C
              CALL PLCHHQ (CFUX(.97),CFUY(.950),'ORIGINAL RECTANGULAR GR
     +ID',.012,0.,1.)
C
              CALL PLCHHQ (CFUX(.97),CFUY(.928),
     +                     'Portions of the grid with missing data',
     +                                                       .010,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.908),
     +                     '(over land areas) are drawn',
     +                                                       .010,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.888),
     +                     'in a lighter gray.',
     +                                                       .010,0.,1.)
C
              CALL PLCHHQ (CFUX(.03),CFUY(.084),'Grid is gray.',
     +                                                      .012,0.,-1.)
              CALL PLCHHQ (CFUX(.03),CFUY(.060),'Shorelines are blue.',
     +                                                      .012,0.,-1.)
              CALL PLCHHQ (CFUX(.03),CFUY(.036),'Parallels/meridians are
     + red.',.012,0.,-1.)
C
C Advance the frame.
C
              CALL FRAME
C
            END IF
C
C Triangular mesh (with edges between blocked triangles somewhat fainter
C than the rest):
C
            IF (MOD(IMSH/4,2).NE.0) THEN
C
              PRINT * , 'DRAWING TRIANGULAR MESH'
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 MAPLOT
C
C Label the third frame.
C
              CALL PLCHHQ (CFUX(.03),CFUY(.946),'ORCA GRID',.024,0.,-1.)
C
              IF (     IDIR.EQ.1) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 1',
     +                                                      .012,0.,-1.)
              ELSE IF (IDIR.EQ.2) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 2',
     +                                                      .012,0.,-1.)
              ELSE IF (IDIR.EQ.3) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 3',
     +                                                      .012,0.,-1.)
              ELSE IF (IDIR.EQ.4) THEN
                CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 4',
     +                                                      .012,0.,-1.)
              END IF
C
              CALL PLCHHQ (CFUX(.97),CFUY(.950),'DERIVED TRIANGULAR MESH
     +',.012,0.,1.)
C
              CALL PLCHHQ (CFUX(.97),CFUY(.928),
     +                     'Triangles from rectangles over',
     +                                                       .010,0.,1.)
              CALL PLCHHQ (CFUX(.97),CFUY(.908),
     +                     'land have been omitted.',
     +                                                       .010,0.,1.)
C
              CALL PLCHHQ (CFUX(.03),CFUY(.084),'Mesh is gray.',
     +                                                      .012,0.,-1.)
              CALL PLCHHQ (CFUX(.03),CFUY(.060),'Shorelines are blue.',
     +                                                      .012,0.,-1.)
              CALL PLCHHQ (CFUX(.03),CFUY(.036),'Parallels/meridians are
     + red.',.012,0.,-1.)
C
C Advance the frame.
C
              CALL FRAME
C
            END IF
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 MAPLOT
C
C Label the fourth frame.
C
            CALL PLCHHQ (CFUX(.03),CFUY(.946),'ORCA GRID',.024,0.,-1.)
C
            IF (     IDIR.EQ.1) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 1',
     +                                                      .012,0.,-1.)
            ELSE IF (IDIR.EQ.2) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 2',
     +                                                      .012,0.,-1.)
            ELSE IF (IDIR.EQ.3) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 3',
     +                                                      .012,0.,-1.)
            ELSE IF (IDIR.EQ.4) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 4',
     +                                                      .012,0.,-1.)
            END IF
C
            CALL PLCHHQ (CFUX(.97),CFUY(.950),'SIMPLE CONTOURS ON',
     +                                                       .012,0.,1.)
            CALL PLCHHQ (CFUX(.97),CFUY(.926),'TRIANGULAR MESH',
     +                                                       .012,0.,1.)
            CALL PLCHHQ (CFUX(.97),CFUY(.904),
     +                   'Actual temperature',
     +                                                       .010,0.,1.)
            CALL PLCHHQ (CFUX(.97),CFUY(.884),
     +                   'data are used.',
     +                                                       .010,0.,1.)
C
            CALL PLCHHQ (CFUX(.03),CFUY(.060),'Shorelines are blue.',
     +                                                      .012,0.,-1.)
            CALL PLCHHQ (CFUX(.03),CFUY(.036),'Parallels/meridians are r
     +ed.',.012,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 MAPLOT
C
C Label the fifth frame.
C
            CALL PLCHHQ (CFUX(.03),CFUY(.946),'ORCA GRID',
     +                                                      .024,0.,-1.)
C
            IF (     IDIR.EQ.1) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 1',
     +                                                      .012,0.,-1.)
            ELSE IF (IDIR.EQ.2) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 2',
     +                                                      .012,0.,-1.)
            ELSE IF (IDIR.EQ.3) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 3',
     +                                                      .012,0.,-1.)
            ELSE IF (IDIR.EQ.4) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 4',
     +                                                      .012,0.,-1.)
            END IF
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.)
            CALL PLCHHQ (CFUX(.97),CFUY(.904),
     +                   'Actual temperature',
     +                                                       .010,0.,1.)
            CALL PLCHHQ (CFUX(.97),CFUY(.884),
     +                   'data are used.',
     +                                                       .010,0.,1.)
C
            CALL PLCHHQ (CFUX(.03),CFUY(.132),'Off-mesh',
     +                                                      .012,0.,-1.)
            CALL PLCHHQ (CFUX(.03),CFUY(.108),'areas of the',
     +                                                      .012,0.,-1.)
            CALL PLCHHQ (CFUX(.03),CFUY(.084),'globe are yellow.',
     +                                                      .012,0.,-1.)
            CALL PLCHHQ (CFUX(.03),CFUY(.060),'Shorelines are blue.',
     +                                                      .012,0.,-1.)
            CALL PLCHHQ (CFUX(.03),CFUY(.036),'Parallels/meridians are r
     +ed.',.012,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'
            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 PLOTIT (0,0,2)
            CALL GSPLCI (13)
            CALL MAPGRD
            CALL PLOTIT (0,0,2)
            CALL GSPLCI (14)
            CALL MAPLOT
C
C Label the sixth frame.
C
            CALL PLCHHQ (CFUX(.03),CFUY(.946),'ORCA GRID',.024,0.,-1.)
C
            IF (     IDIR.EQ.1) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 1',
     +                                                      .012,0.,-1.)
            ELSE IF (IDIR.EQ.2) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 2',
     +                                                      .012,0.,-1.)
            ELSE IF (IDIR.EQ.3) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 3',
     +                                                      .012,0.,-1.)
            ELSE IF (IDIR.EQ.4) THEN
              CALL PLCHHQ (CFUX(.03),CFUY(.874),'VIEWPOINT 4',
     +                                                      .012,0.,-1.)
            END IF
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),
     +                   'Actual temperature',
     +                                                       .010,0.,1.)
            CALL PLCHHQ (CFUX(.97),CFUY(.884),
     +                   'data are used.',
     +                                                       .010,0.,1.)
C
            CALL PLCHHQ (CFUX(.03),CFUY(.132),'Off-mesh',
     +                                                      .012,0.,-1.)
            CALL PLCHHQ (CFUX(.03),CFUY(.108),'areas of the',
     +                                                      .012,0.,-1.)
            CALL PLCHHQ (CFUX(.03),CFUY(.084),'globe are yellow.',
     +                                                      .012,0.,-1.)
            CALL PLCHHQ (CFUX(.03),CFUY(.060),'Shorelines are blue.',
     +                                                      .012,0.,-1.)
            CALL PLCHHQ (CFUX(.03),CFUY(.036),'Parallels/meridians are r
     +ed.',.012,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 GTORCA (RPNT,MPNT,NPNT,LOPN,
     +                   IEDG,MEDG,NEDG,LOEN,
     +                   ITRI,MTRI,NTRI,LOTN,
     +                             CLAT,CLON)
C
        DIMENSION RPNT(MPNT),IEDG(MEDG),ITRI(MTRI)
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 The ORCA grid is rectangular; at each point, we have a latitude, a
C longitude, and a data value.  From this, we construct a triangular
C mesh.
C
C Define some dimensions.
C
        PARAMETER (IDIM=181,JDIM=148)
        PARAMETER (IDM1=IDIM-1,JDM1=JDIM-1)
        PARAMETER (IDP1=IDIM+1,JDP1=JDIM+1)
C
C Dimension arrays to receive the original data from the NetCDF file.
C RLAT is an array of latitudes, RLON an array of longitudes, and
C DVAL an array of data.  Each of these arrays has two points of
C overlap in the first dimension: for each value of J from 1 to JDP1,
C the grid points indexed [1,J] and [2,J] are the same as the ones
C indexed [IDIM,J] and [IDP1,J], respectively.  There are also some
C points of overlap in the second dimension: the grid points with
C indices [I,JDP1], for I from 1 to IDIM, trace out a long narrow
C ellipse near the North Pole which is a duplicate of one traced out
C by the grid points with indices [I,JDM1], for I from 1 to IDIM, but
C in reversed order; the grid points with indices [I,JDIM], for I from
C 1 to IDIM, trace out a degenerate ellipse with a minor axis of zero
C (a line traced from beginning to end and back again).
C
        DIMENSION RLAT(IDP1,JDP1),RLON(IDP1,JDP1),DVAL(IDP1,JDP1,5)
C
C Declare arrays to receive a subset of the ORCA grid data.  XLAT is the
C array of latitudes, XLON is the array of longitudes, and ZDAT is an
C array of data.  One "horizontal" line (for J=IDP1) and one "vertical"
C line (for I=IDP1) of the original ORCA grid data are discarded in
C the copy to these arrays, but each of them still has one point of
C overlap in the first dimension: for each value of J from 1 to JDIM,
C the grid point indexed [1,J] is the same as the one indexed [IDIM,J].
C
        COMMON /CMORCA/ XLAT(IDIM,JDIM),XLON(IDIM,JDIM),ZDAT(IDIM,JDIM)
        SAVE   /CMORCA/
C
C Each cell of the rectangular grid for which data is available at all
C four corners is split in half along one of its diagonals; resulting
C triangles form the triangular mesh.  As we create the triangular mesh,
C we have to keep track of where the points and edges of the rectangular
C grid were put (so as to avoid the problem of duplicating points and
C edges in the structure).  The array ISCR is used for this.
C
        DIMENSION ISCR(4,IDIM,JDIM)
C
C Declare external a routine to tell CTTMRG about points of overlap on
C the ORCA grid.
C
        EXTERNAL MIORCA
C
C Define the value that indicates a missing datum in the ORCA grid.
C
        DATA SVAL / 1.E20 /
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('ctorca.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.
C
C       ISTA=NF_INQ_VARID(NCID,'nav_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,RLAT)
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.
C
C       ISTA=NF_INQ_VARID(NCID,'nav_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,RLON)
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 data.
C
C       ISTA=NF_INQ_VARID(NCID,'sosstsst',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,DVAL)
C       IF (ISTA.NE.NF_NOERR) THEN
C         PRINT * , 'ERROR RETURN FROM NF_GET_VAR_REAL: ',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 "ctorca.dat".)
C
C       WRITE (11,'(5E16.8)') RLAT
C       WRITE (11,'(5E16.8)') RLON
C       WRITE (11,'(5E16.8)') DVAL
C
C Read the required data from the ASCII file "ctorca.dat".
C
        OPEN (11,FILE='ctorca.dat',STATUS='OLD',FORM='FORMATTED')
C
        READ (11,'(5E16.8)') RLAT
        READ (11,'(5E16.8)') RLON
        READ (11,'(5E16.8)') DVAL
C
        CLOSE (11)
C
C Copy the data from the original arrays to the subset arrays.
C
        DO 102 I=1,IDIM
          DO 101 J=1,JDIM
            XLAT(I,J)=RLAT(I+1,J)
            XLON(I,J)=RLON(I+1,J)
            ZDAT(I,J)=DVAL(I+1,J,1)
  101     CONTINUE
  102   CONTINUE
C
C Call a general-purpose subroutine that accepts a rectangular grid
C mapped onto the surface of the globe and returns a triangular mesh
C equivalent to it.
C
        CALL CTTMRG (IDIM,JDIM,XLAT,XLON,ZDAT,ISCR,SVAL,MIORCA,
     +               RPNT,MPNT,NPNT,LOPN,
     +               IEDG,MEDG,NEDG,LOEN,
     +               ITRI,MTRI,NTRI,LOTN)
C
C Return the latitude and longitude of the approximate center point of
C the mesh on the globe.
C
        CLAT=XLAT(IDIM/2,JDIM/2)
        CLON=XLON(IDIM/2,JDIM/2)
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE MIORCA (IDIM,JDIM,IINI,JINI,IINO,JINO)
C
C Given the dimensions, IDIM and JDIM, of a lat/lon grid on the globe,
C and the indices, IINI and JINI, of a point on the grid, this routine
C returns the indices, IINO and JINO, of that coincident point on the
C grid which is to be used to represent it.  This version assumes that
C the first and last rows of the grid each map into a single point
C (perhaps the south pole and the north pole, respectively) and that
C the right and left edges of the grid lie on top of each other.
C
        IF      (JINI.EQ.JDIM) THEN  !  top row of grid
          IINO=MIN(IINI,IDIM+1-IINI)
          JINO=JDIM
        ELSE IF (IINI.EQ.IDIM) THEN  !  point in last column of grid
          IINO=1
          JINO=JINI
        ELSE                         !  all other points of the grid
          IINO=IINI
          JINO=JINI
        END IF
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE DSORCA
C
C Draw the original ORCA grid, using data passed in a common block.
C
        PARAMETER (IDIM=181,JDIM=148,IDM1=IDIM-1,JDM1=JDIM-1)
C
        COMMON /CMORCA/ XLAT(IDIM,JDIM),XLON(IDIM,JDIM),ZDAT(IDIM,JDIM)
        SAVE   /CMORCA/
C
C Define the value that indicates a missing datum in the ORCA grid.
C
        DATA SVAL / 1.E20 /
C
C Draw the grid box by box, skipping any box with a side longer than
C 2.5 degrees and any box in the list embedded in the function KILLIT.
C
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (8)
C
        DO 102 I=1,IDM1
          DO 101 J=1,JDM1
            IF (ADSGCR(XLAT(I  ,J  ),XLON(I  ,J  ),
     +                 XLAT(I+1,J  ),XLON(I+1,J  )).GT.2.5) GO TO 101
            IF (ADSGCR(XLAT(I+1,J  ),XLON(I+1,J  ),
     +                 XLAT(I+1,J+1),XLON(I+1,J+1)).GT.2.5) GO TO 101
            IF (ADSGCR(XLAT(I+1,J+1),XLON(I+1,J+1),
     +                 XLAT(I  ,J+1),XLON(I  ,J+1)).GT.2.5) GO TO 101
            IF (ADSGCR(XLAT(I  ,J+1),XLON(I  ,J+1),
     +                 XLAT(I  ,J  ),XLON(I  ,J  )).GT.2.5) GO TO 101
            IF (KILLIT(I,J).NE.0) GO TO 101
            IF (ZDAT(I  ,J  ).NE.SVAL.AND.ZDAT(I  ,J+1).NE.SVAL) THEN
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (8)
            ELSE
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (9)
            END IF
            CALL DRSGCR (XLAT(I  ,J  ),XLON(I  ,J  ),
     +                   XLAT(I  ,J+1),XLON(I  ,J+1))
            IF (ZDAT(I  ,J+1).NE.SVAL.AND.ZDAT(I+1,J+1).NE.SVAL) THEN
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (8)
            ELSE
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (9)
            END IF
            CALL DRSGCR (XLAT(I  ,J+1),XLON(I  ,J+1),
     +                   XLAT(I+1,J+1),XLON(I+1,J+1))
            IF (ZDAT(I+1,J+1).NE.SVAL.AND.ZDAT(I+1,J  ).NE.SVAL) THEN
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (8)
            ELSE
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (9)
            END IF
            CALL DRSGCR (XLAT(I+1,J+1),XLON(I+1,J+1),
     +                   XLAT(I+1,J  ),XLON(I+1,J  ))
            IF (ZDAT(I+1,J  ).NE.SVAL.AND.ZDAT(I  ,J  ).NE.SVAL) THEN
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (8)
            ELSE
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (9)
            END IF
            CALL DRSGCR (XLAT(I+1,J  ),XLON(I+1,J  ),
     +                   XLAT(I  ,J  ),XLON(I  ,J  ))
  101     CONTINUE
  102   CONTINUE
C
C Draw the continental outlines and lat/lon grid.
C
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (13)
        CALL MAPGRD
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (14)
        CALL MAPLOT
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE DCORCA
C
C Draw the original ORCA grid, using data passed in a common block.
C
        PARAMETER (IDIM=181,JDIM=148,IDM1=IDIM-1,JDM1=JDIM-1)
C
        COMMON /CMORCA/ XLAT(IDIM,JDIM),XLON(IDIM,JDIM),ZDAT(IDIM,JDIM)
        SAVE   /CMORCA/
C
        COMMON /CMGRID/ IWID,IFCI,MCOL,NCOL
        SAVE   /CMGRID/
C
C Declare arrays to use for solid fill.
C
        DIMENSION XBOX(4),YBOX(4)
C
C Define the major and minor axes of the ellipses that mark selected
C points on the grid.
C
        DATA RADA,RADB / .010,.010 /
C
C Define a constant to convert from degrees to radians.
C
        DATA DTOR / .017453292519943 /
C
C Fill the grid box by box, skipping any box with a side longer than
C 2.5 degrees and any box in the list embedded in the function KILLIT.
C
        DO 102 I=1,IDM1
          ICOL=INT(10.*REAL(I-1)/REAL(IDM1))
          DO 101 J=1,JDM1
            JCOL=INT(10.*REAL(J-1)/REAL(JDM1))
            CALL PLOTIT (0,0,2)
            CALL GSPLCI (51+10*ICOL+JCOL)
            CALL GSFACI (51+10*ICOL+JCOL)
            IF (ADSGCR(XLAT(I  ,J  ),XLON(I  ,J  ),
     +                 XLAT(I+1,J  ),XLON(I+1,J  )).GT.2.5) GO TO 101
            IF (ADSGCR(XLAT(I+1,J  ),XLON(I+1,J  ),
     +                 XLAT(I+1,J+1),XLON(I+1,J+1)).GT.2.5) GO TO 101
            IF (ADSGCR(XLAT(I+1,J+1),XLON(I+1,J+1),
     +                 XLAT(I  ,J+1),XLON(I  ,J+1)).GT.2.5) GO TO 101
            IF (ADSGCR(XLAT(I  ,J+1),XLON(I  ,J+1),
     +                 XLAT(I  ,J  ),XLON(I  ,J  )).GT.2.5) GO TO 101
            IF (KILLIT(I,J).NE.0) GO TO 101
            XCOA=COS(DTOR*XLAT(I  ,J  ))*COS(DTOR*XLON(I  ,J  ))
            YCOA=COS(DTOR*XLAT(I  ,J  ))*SIN(DTOR*XLON(I  ,J  ))
            ZCOA=SIN(DTOR*XLAT(I  ,J  ))
            CALL CTMXYZ (1,XCOA,YCOA,ZCOA,XBOX(1),YBOX(1))
            IF (XBOX(1).EQ.1.E12.OR.YBOX(1).EQ.1.E12) GO TO 101
            XCOA=COS(DTOR*XLAT(I+1,J  ))*COS(DTOR*XLON(I+1,J  ))
            YCOA=COS(DTOR*XLAT(I+1,J  ))*SIN(DTOR*XLON(I+1,J  ))
            ZCOA=SIN(DTOR*XLAT(I+1,J  ))
            CALL CTMXYZ (1,XCOA,YCOA,ZCOA,XBOX(2),YBOX(2))
            IF (XBOX(2).EQ.1.E12.OR.YBOX(2).EQ.1.E12) GO TO 101
            XCOA=COS(DTOR*XLAT(I+1,J+1))*COS(DTOR*XLON(I+1,J+1))
            YCOA=COS(DTOR*XLAT(I+1,J+1))*SIN(DTOR*XLON(I+1,J+1))
            ZCOA=SIN(DTOR*XLAT(I+1,J+1))
            CALL CTMXYZ (1,XCOA,YCOA,ZCOA,XBOX(3),YBOX(3))
            IF (XBOX(3).EQ.1.E12.OR.YBOX(3).EQ.1.E12) GO TO 101
            XCOA=COS(DTOR*XLAT(I  ,J+1))*COS(DTOR*XLON(I  ,J+1))
            YCOA=COS(DTOR*XLAT(I  ,J+1))*SIN(DTOR*XLON(I  ,J+1))
            ZCOA=SIN(DTOR*XLAT(I  ,J+1))
            CALL CTMXYZ (1,XCOA,YCOA,ZCOA,XBOX(4),YBOX(4))
            IF (XBOX(4).EQ.1.E12.OR.YBOX(4).EQ.1.E12) GO TO 101
            CALL GFA (4,XBOX,YBOX)
  101     CONTINUE
  102   CONTINUE
C
C Draw the continental outlines and lat/lon grid.
C
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (13)
        CALL MAPGRD
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (14)
        CALL MAPLOT
C
C In gray, draw those interior lines of the grid that separate boxes of
C different colors.
C
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (11)
C
        DO 104 I=2,IDM1
          DO 103 J=1,JDM1
            IF (ADSGCR(XLAT(I,J  ),XLON(I,J  ),
     +                 XLAT(I,J+1),XLON(I,J+1)).GT.2.5) GO TO 103
            IF (INT(10.*REAL(I-1)/REAL(IDM1)).EQ.
     +          INT(10.*REAL(I-2)/REAL(IDM1))) GO TO 103
            CALL DRSGCR (XLAT(I,J),XLON(I,J),XLAT(I,J+1),XLON(I,J+1))
  103     CONTINUE
  104   CONTINUE
C
        DO 106 J=2,JDM1
          DO 105 I=1,IDM1
            IF (ADSGCR(XLAT(I  ,J),XLON(I  ,J),
     +                 XLAT(I+1,J),XLON(I+1,J)).GT.2.5) GO TO 105
            IF (INT(10.*REAL(J-1)/REAL(JDM1)).EQ.
     +          INT(10.*REAL(J-2)/REAL(JDM1))) GO TO 105
            CALL DRSGCR (XLAT(I,J),XLON(I,J),XLAT(I+1,J),XLON(I+1,J))
  105     CONTINUE
  106   CONTINUE
C
C In the foreground color, outline the grid (using a thickened line).
C
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (1)
        CALL GSLWSC (2.)
C
        DO 108 I=1,IDIM,IDIM-1
          DO 107 J=1,JDM1
            IF (ADSGCR(XLAT(I,J  ),XLON(I,J  ),
     +                 XLAT(I,J+1),XLON(I,J+1)).GT.2.5) GO TO 107
            CALL DRSGCR (XLAT(I,J),XLON(I,J),XLAT(I,J+1),XLON(I,J+1))
  107     CONTINUE
  108   CONTINUE
C
        DO 110 J=1,JDIM,JDIM-1
          DO 109 I=1,IDM1
            IF (ADSGCR(XLAT(I  ,J),XLON(I  ,J),
     +                 XLAT(I+1,J),XLON(I+1,J)).GT.2.5) GO TO 109
            CALL DRSGCR (XLAT(I,J),XLON(I,J),XLAT(I+1,J),XLON(I+1,J))
  109     CONTINUE
  110   CONTINUE
C
        CALL PLOTIT (0,0,2)
        CALL GSLWSC (1.)
C
C Mark nine reference points on the grid.
C
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (1)
        CALL GSFACI (1)
C
        CALL MAPTRA (XLAT(1,JDIM),XLON(1,JDIM),XVAL,YVAL)
        IF (XVAL.NE.1.E12.AND.YVAL.NE.1.E12) THEN
          CALL ELLIPS (CUFX(XVAL),CUFY(YVAL),RADA,RADB,0.,4.)
        END IF
C
        CALL MAPTRA (XLAT(IDIM,JDIM),XLON(IDIM,JDIM),XVAL,YVAL)
        IF (XVAL.NE.1.E12.AND.YVAL.NE.1.E12) THEN
          CALL ELLIPS (CUFX(XVAL),CUFY(YVAL),RADA,RADB,0.,4.)
        END IF
C
        CALL MAPTRA (XLAT(1,1),XLON(1,1),XVAL,YVAL)
        IF (XVAL.NE.1.E12.AND.YVAL.NE.1.E12) THEN
          CALL ELLIPS (CUFX(XVAL),CUFY(YVAL),RADA,RADB,0.,4.)
        END IF
C
        CALL MAPTRA (XLAT(IDIM,1),XLON(IDIM,1),XVAL,YVAL)
        IF (XVAL.NE.1.E12.AND.YVAL.NE.1.E12) THEN
          CALL ELLIPS (CUFX(XVAL),CUFY(YVAL),RADA,RADB,0.,4.)
        END IF
C
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (1)
        CALL GSFACI (0)
C
        IF (MOD(IDIM,2).EQ.0) THEN
          CALL FPSGCR (XLAT(IDIM/2  ,JDIM),XLON(IDIM/2  ,JDIM),
     +                 XLAT(IDIM/2+1,JDIM),XLON(IDIM/2+1,JDIM),
     +                                                     .5,CLAT,CLON)
        ELSE
          CLAT=XLAT(IDIM/2+1,JDIM)
          CLON=XLON(IDIM/2+1,JDIM)
        END IF
        CALL MAPTRA (CLAT,CLON,XVAL,YVAL)
        IF (XVAL.NE.1.E12.AND.YVAL.NE.1.E12) THEN
          CALL ELLIPS (CUFX(XVAL),CUFY(YVAL),RADA,RADB,0.,4.)
        END IF
C
        IF (MOD(JDIM,2).EQ.0) THEN
          CALL FPSGCR (XLAT(1,JDIM/2  ),XLON(1,JDIM/2  ),
     +                 XLAT(1,JDIM/2+1),XLON(1,JDIM/2+1),
     +                                                     .5,CLAT,CLON)
        ELSE
          CLAT=XLAT(1,JDIM/2+1)
          CLON=XLON(1,JDIM/2+1)
        END IF
        CALL MAPTRA (CLAT,CLON,XVAL,YVAL)
        IF (XVAL.NE.1.E12.AND.YVAL.NE.1.E12) THEN
          CALL ELLIPS (CUFX(XVAL),CUFY(YVAL),RADA,RADB,0.,4.)
        END IF
C
        IF (MOD(JDIM,2).EQ.0) THEN
          CALL FPSGCR (XLAT(IDIM,JDIM/2  ),XLON(IDIM,JDIM/2  ),
     +                 XLAT(IDIM,JDIM/2+1),XLON(IDIM,JDIM/2+1),
     +                                                     .5,CLAT,CLON)
        ELSE
          CLAT=XLAT(IDIM,JDIM/2+1)
          CLON=XLON(IDIM,JDIM/2+1)
        END IF
        CALL MAPTRA (CLAT,CLON,XVAL,YVAL)
        IF (XVAL.NE.1.E12.AND.YVAL.NE.1.E12) THEN
          CALL ELLIPS (CUFX(XVAL),CUFY(YVAL),RADA,RADB,0.,4.)
        END IF
C
        IF (MOD(IDIM,2).EQ.0) THEN
          CALL FPSGCR (XLAT(IDIM/2  ,1),XLON(IDIM/2  ,1),
     +                 XLAT(IDIM/2+1,1),XLON(IDIM/2+1,1),
     +                                                     .5,CLAT,CLON)
        ELSE
          CLAT=XLAT(IDIM/2+1,1)
          CLON=XLON(IDIM/2+1,1)
        END IF
        CALL MAPTRA (CLAT,CLON,XVAL,YVAL)
        IF (XVAL.NE.1.E12.AND.YVAL.NE.1.E12) THEN
          CALL ELLIPS (CUFX(XVAL),CUFY(YVAL),RADA,RADB,0.,4.)
        END IF
C
        CALL PLOTIT (0,0,2)
        CALL GSPLCI (1)
        CALL GSFACI (4)
C
        IF (MOD(IDIM,2).EQ.0) THEN
          IF (MOD(JDIM,2).EQ.0) THEN
            CALL IPIQOS(XLAT(IDIM/2  ,JDIM/2  ),XLON(IDIM/2  ,JDIM/2  ),
     +                  XLAT(IDIM/2+1,JDIM/2  ),XLON(IDIM/2+1,JDIM/2  ),
     +                  XLAT(IDIM/2  ,JDIM/2+1),XLON(IDIM/2  ,JDIM/2+1),
     +                  XLAT(IDIM/2+1,JDIM/2+1),XLON(IDIM/2+1,JDIM/2+1),
     +                                                  .5,.5,CLAT,CLON)
          ELSE
            CALL FPSGCR(XLAT(IDIM/2  ,JDIM/2+1),XLON(IDIM/2  ,JDIM/2+1),
     +                  XLAT(IDIM/2+1,JDIM/2+1),XLON(IDIM/2+1,JDIM/2+1),
     +                                                     .5,CLAT,CLON)
          END IF
        ELSE
          IF (MOD(JDIM,2).EQ.0) THEN
            CALL FPSGCR(XLAT(IDIM/2+1,JDIM/2  ),XLON(IDIM/2+1,JDIM/2  ),
     +                  XLAT(IDIM/2+1,JDIM/2+1),XLON(IDIM/2+1,JDIM/2+1),
     +                                                    .5,CLAT,CLON)
          ELSE
           CLAT=XLAT(IDIM/2+1,JDIM/2+1)
           CLON=XLON(IDIM/2+1,JDIM/2+1)
          END IF
        END IF
        CALL MAPTRA (CLAT,CLON,XVAL,YVAL)
        IF (XVAL.NE.1.E12.AND.YVAL.NE.1.E12) THEN
          CALL ELLIPS (CUFX(XVAL),CUFY(YVAL),RADA,RADB,0.,4.)
        END IF
C
C Draw a shaded reference grid in the lower left corner of the frame.
C
        CALL SHDGRD (.03,.33,.03,.27,IDM1,JDM1,3,7)
C
C Done.
C
        RETURN
C
      END


      FUNCTION KILLIT (I,J)
C
C This function has a non-zero value if and only if the ORCA grid cell
C whose lower left corner has indices I and J should be deleted.
C
        IF ( (I.GE.  7.AND.I.LE.  7.AND.J.GE. 96.AND.J.LE.112).OR.
     +       (I.GE.139.AND.I.LE.139.AND.J.GE.103.AND.J.LE.103).OR.
     +       (I.GE.138.AND.I.LE.139.AND.J.GE. 98.AND.J.LE. 99).OR.
     +       (I.GE.139.AND.I.LE.139.AND.J.GE.104.AND.J.LE.104).OR.
     +       (I.GE.141.AND.I.LE.142.AND.J.GE.106.AND.J.LE.111).OR.
     +       (I.GE.143.AND.I.LE.144.AND.J.GE.109.AND.J.LE.110).OR.
     +       (I.GE.143.AND.I.LE.143.AND.J.GE.105.AND.J.LE.105).OR.
     +       (I.GE.163.AND.I.LE.165.AND.J.GE. 85.AND.J.LE. 85)) THEN
          KILLIT=1
        ELSE
          KILLIT=0
        END IF
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 DCFGRD (IAR1,IAR2,IAR3,IAR4,CCLL,CCLR,CCUL,CCUR)
C
        DIMENSION CCLL(3),CCLR(3),CCUL(3),CCUR(3),CCGC(3)
C
C This routine defines the colors to be used by SHDGRD.  CCLL, CCLR,
C CCUL, and CCUR define the RGB color components for the lower-left
C corner, the lower-right corner, the upper-left corner, and the
C upper-right corner, respectively, of the grid.  Other colors to be
C used are defined by interpolation.
C
        COMMON /CMGRID/ IWID,IFCI,MCOL,NCOL
        SAVE   /CMGRID/
C
C Transfer first four arguments to common block shared with SHDGRD.
C
        IWID=IAR1
        IFCI=IAR2
        MCOL=IAR3
        NCOL=IAR4
C
        DO 103 I=0,MCOL-1
          P=REAL(I)/REAL(MCOL-1)
          DO 102 J=0,NCOL-1
            Q=REAL(J)/REAL(NCOL-1)
            DO 101 K=1,3
              CCGC(K)=(1.-P)*((1.-Q)*CCLL(K)+Q*CCUL(K))+
     +                    P *((1.-Q)*CCLR(K)+Q*CCUR(K))
              CALL GSCR (IWID,IFCI+NCOL*I+J,CCGC(1),CCGC(2),CCGC(3))
  101       CONTINUE
  102     CONTINUE
  103   CONTINUE
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE SHDGRD (XVPL,XVPR,YVPB,YVPT,IMAX,JMAX,ITYP,IMRK)
C
C This routine draws a representation of a grid of IMAX by JMAX boxes
C (having IMAX+1 vertical lines and JMAX+1 horizontal lines) in the part
C of the plotter frame specified by XVPL, XVPR, YVPB, and YVPT.  Colors
C used are those specified by a superimposed grid of MCOL by NCOL colors
C with indices IFCI through IFCI+MCOL*NCOL-1.  If ITYP is 0, the grid
C is drawn, but not filled.  If ITYP is 1, the grid is filled.  If ITYP
C is 2, the grid is filled and then outlined in the foreground color.
C If ITYP is 3, the grid is filled, some of the interior grid lines (the
C ones between boxes of different colors) are drawn in gray, and then
C the grid is outlined in the foreground color.  If ITYP is 4, the grid
C is filled, all of the interior grid lines are drawn in gray, and then
C the grid is outlined in the foreground color.  If IMRK is 0, no points
C of the grid are marked with filled circles; if IMRK is 1, the corner
C points are marked with black-filled circles; if IMRK is 2, the edge
C points are marked with white-filled black circles; if IMRK is 4, the
C center point is marked with a red-filled black circle; sums of these
C values can be used to get combinations of behaviors.
C
        COMMON /CMGRID/ IWID,IFCI,MCOL,NCOL
        SAVE   /CMGRID/
C
        DIMENSION XBOX(5),YBOX(5)
C
C Define the major and minor axes of the ellipses that mark selected
C points on the grid.
C
        DATA RADA,RADB / .005,.005 /
C
C Save the current SET call.
C
        CALL GETSET (EVPL,EVPR,EVPB,EVPT,EWDL,EWDR,EWDB,EWDT,LNLG)
C
C Compute the coordinates of the center of the viewport in which the
C grid is to be drawn, its half-height, and its half-width.
C
        XVPC=(XVPR+XVPL)/2.
        YVPC=(YVPB+YVPT)/2.
        HVPW=(XVPR-XVPL)/2.
        HVPH=(YVPT-YVPB)/2.
C
C Compute values defining the exact viewport in which the grid is to
C be drawn and redo the set call appropriately.
C
        IF (REAL(IMAX)/REAL(JMAX).LE.(XVPR-XVPL)/(YVPT-YVPB)) THEN
          SVPL=XVPC-HVPH*(REAL(IMAX)/REAL(JMAX))
          SVPR=XVPC+HVPH*(REAL(IMAX)/REAL(JMAX))
          SVPB=YVPB
          SVPT=YVPT
        ELSE
          SVPL=XVPL
          SVPR=XVPR
          SVPB=YVPC-HVPW*(REAL(JMAX)/REAL(IMAX))
          SVPT=YVPC+HVPW*(REAL(JMAX)/REAL(IMAX))
        END IF
C
        CALL SET (SVPL,SVPR,SVPB,SVPT,0.,REAL(IMAX),0.,REAL(JMAX),1)
C
C If the lines are to be drawn in different colors ...
C
        IF (ITYP.EQ.0) THEN
C
C ... do it.
C
          DO 102 I=0,IMAX-1
            ICOL=INT(REAL(MCOL)*REAL(I)/REAL(IMAX))
            DO 101 J=0,JMAX-1
              JCOL=INT(REAL(NCOL)*REAL(J)/REAL(JMAX))
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (IFCI+NCOL*ICOL+JCOL)
              XBOX(1)=REAL(I)
              XBOX(2)=REAL(I+1)
              XBOX(3)=REAL(I+1)
              XBOX(4)=REAL(I)
              XBOX(5)=REAL(I)
              YBOX(1)=REAL(J)
              YBOX(2)=REAL(J)
              YBOX(3)=REAL(J+1)
              YBOX(4)=REAL(J+1)
              YBOX(5)=REAL(J)
              CALL GPL (5,XBOX,YBOX)
  101       CONTINUE
  102     CONTINUE
C
C Otherwise, the boxes are to be filled in different colors, so ...
C
        ELSE
C
C ... do it.
C
          DO 104 I=0,IMAX-1
            ICOL=INT(REAL(MCOL)*REAL(I)/REAL(IMAX))
            DO 103 J=0,JMAX-1
              JCOL=INT(REAL(NCOL)*REAL(J)/REAL(JMAX))
              CALL GSFACI (IFCI+NCOL*ICOL+JCOL)
              XBOX(1)=REAL(I)
              XBOX(2)=REAL(I+1)
              XBOX(3)=REAL(I+1)
              XBOX(4)=REAL(I)
              YBOX(1)=REAL(J)
              YBOX(2)=REAL(J)
              YBOX(3)=REAL(J+1)
              YBOX(4)=REAL(J+1)
              CALL GFA (4,XBOX,YBOX)
  103       CONTINUE
  104     CONTINUE
C
C In addition, if specified, draw lines on the grid ...
C
          IF (ITYP.GT.1) THEN
C
C ... either just the border ...
C
            IF (ITYP.GT.2) THEN
C
C ... or all the lines ...
C
              CALL PLOTIT (0,0,2)
              CALL GSPLCI (11)
C
              DO 105 I=1,IMAX-1
                IF (ITYP.GT.3.OR.
     +              INT(REAL(MCOL)*REAL(I-1)/REAL(IMAX)).NE.
     +              INT(REAL(MCOL)*REAL(I  )/REAL(IMAX)))
     +            CALL LINE (REAL(I),0.,REAL(I),REAL(JMAX))
  105         CONTINUE
C
              DO 106 J=1,JMAX-1
                IF (ITYP.GT.3.OR.
     +              INT(REAL(NCOL)*REAL(J-1)/REAL(JMAX)).NE.
     +              INT(REAL(NCOL)*REAL(J  )/REAL(JMAX)))
     +          CALL LINE (0.,REAL(J),REAL(IMAX),REAL(J))
  106         CONTINUE
C
            END IF
C
            CALL PLOTIT (0,0,2)
            CALL GSPLCI (1)
            CALL GSLWSC (2.)
C
            CALL LINE (        0.,        0.,REAL(IMAX),        0.)
            CALL LINE (REAL(IMAX),        0.,REAL(IMAX),REAL(JMAX))
            CALL LINE (REAL(IMAX),REAL(JMAX),        0.,REAL(JMAX))
            CALL LINE (        0.,REAL(JMAX),        0.,        0.)
C
            CALL PLOTIT (0,0,2)
            CALL GSLWSC (1.)
C
          END IF
C
        END IF
C
C If specified, mark a selected subset of nine different points of
C reference on the grid.
C
        IF (MOD(IMRK,2).NE.0) THEN
C
          CALL PLOTIT (0,0,2)
          CALL GSPLCI (1)
          CALL GSFACI (1)
C
          CALL ELLIPS (SVPL,SVPT,RADA,RADB,0.,4.)
          CALL ELLIPS (SVPR,SVPT,RADA,RADB,0.,4.)
          CALL ELLIPS (SVPL,SVPB,RADA,RADB,0.,4.)
          CALL ELLIPS (SVPR,SVPB,RADA,RADB,0.,4.)
C
        END IF
C
        IF (MOD(IMRK/2,2).NE.0) THEN
C
          CALL PLOTIT (0,0,2)
          CALL GSPLCI (1)
          CALL GSFACI (0)
C
          CALL ELLIPS (.5*(SVPL+SVPR),SVPT,RADA,RADB,0.,4.)
          CALL ELLIPS (SVPL,.5*(SVPB+SVPT),RADA,RADB,0.,4.)
          CALL ELLIPS (SVPR,.5*(SVPB+SVPT),RADA,RADB,0.,4.)
          CALL ELLIPS (.5*(SVPL+SVPR),SVPB,RADA,RADB,0.,4.)
C
        END IF
C
        IF (MOD(IMRK/4,2).NE.0) THEN
C
          CALL PLOTIT (0,0,2)
          CALL GSPLCI (1)
          CALL GSFACI (4)
C
          CALL ELLIPS (.5*(SVPL+SVPR),.5*(SVPB+SVPT),RADA,RADB,0.,4.)
C
        END IF
C
        CALL GSFACI (1)
C
C Label the bottom of the grid.
C
        CALL PLCHHQ (CFUX(.5*(SVPL+SVPR)),CFUY(SVPB-.018),
     +                                                'I-->',.012,0.,0.)
C
C Label the left edge of the grid.
C
        CALL PLCHHQ (CFUX(SVPL-.018),CFUY(.5*(SVPB+SVPT)),
     +                                               'J-->',.012,90.,0.)
C
C Restore the original SET call.
C
        CALL SET (EVPL,EVPR,EVPB,EVPT,EWDL,EWDR,EWDB,EWDT,LNLG)
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE ELLIPS (XCFR,YCFR,RADA,RADB,ROTD,DSTP)
C
C This routine fills an ellipse.  The arguments are as follows:
C
C   XCFR and YCFR are the coordinates of the center of the ellipse, in
C   the fractional coordinate system.
C
C   RADA is the length of the semimajor axis of the ellipse (i.e., the
C   distance from the center of the ellipse to one of the two points on
C   the ellipse which are furthest from the center).  This is a distance
C   in the fractional coordinate system.
C
C   RADB is the length of the semiminor axis of the ellipse (i.e., the
C   distance from the center of the ellipse to one of the two points on
C   the ellipse which are nearest to the center).  This is a distance in
C   the fractional coordinate system.
C
C   ROTD is a rotation angle, in degrees.  If ROTD is 0, the major axis
C   of the ellipse is horizontal.  If ROTD is 90, the major axis is
C   vertical.
C
C   DSTP is the step size, in degrees, between any two consecutive
C   points used to draw the ellipse.  The actual value used will be
C   limited to the range from .1 degrees (3600 points used to draw
C   the ellipse) to 90 degrees (4 points used to draw the ellipse).
C
C Declare work arrays to hold the coordinates.
C
        DIMENSION XCRA(3601),YCRA(3601)
C
C DTOR is pi over 180, used to convert an angle from degrees to radians.
C
        DATA DTOR / .017453292519943 /
C
C Get the rotation angle in radians.
C
        ROTR=DTOR*ROTD
C
C Compute the number of steps to be used to draw the ellipse and the
C actual number of degrees for each step.
C
        NSTP=MAX(4,MIN(3600,INT(360./MAX(.1,MIN(90.,DSTP)))))
        RSTP=360./NSTP
C
C Compute coordinates for the ellipse (just some trigonometry).
C
        DO 101 ISTP=0,NSTP
          ANGL=DTOR*REAL(ISTP)*RSTP
          XTMP=RADA*COS(ANGL)
          YTMP=RADB*SIN(ANGL)
          XCRA(ISTP+1)=CFUX(XCFR+XTMP*COS(ROTR)-YTMP*SIN(ROTR))
          YCRA(ISTP+1)=CFUY(YCFR+XTMP*SIN(ROTR)+YTMP*COS(ROTR))
  101   CONTINUE
C
C Fill it.
C
        CALL GFA (NSTP+1,XCRA,YCRA)
C
C Draw it.
C
        CALL GPL (NSTP+1,XCRA,YCRA)
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE IPIQOS (ALAT,ALON,BLAT,BLON,CLAT,CLON,DLAT,DLON,
     +                                       XFRA,YFRA,ELAT,ELON)
C
C (IPIQOS = Interpolate Point In Quadrilateral on Sphere)
C
C This routine, given the latitudes and longitudes of four points (A, B,
C C, and D) forming a "quadrilateral" on the globe and two interpolation
C fractions (XFRA and YFRA, each between 0 and 1, inclusive), finds the
C point E defined by the following diagram and returns its latitude and
C longitude.
C
C                              C------Q----D
C                              |      |    |
C                              |      E    |
C                              |      |    |
C                              |      |    |
C                              A------P----B
C
C P and Q are positioned such that AP/AB = CQ/CD = XFRA and then E is
C positioned such that PE/PQ = YFRA (where "XY" is interpreted to mean
C "the shortest great circle distance from X to Y").
C
C It is assumed that the "quadrilateral" ABDC is "convex" (a working
C definition of which might be that none of the four great circles
C defined by its edges - the ones through A and B, B and D, D and C,
C and C and A - cross it anywhere.  However, this is not verified.
C
C The code is easy:
C
        CALL FPSGCR (ALAT,ALON,BLAT,BLON,XFRA,PLAT,PLON)
        CALL FPSGCR (CLAT,CLON,DLAT,DLON,XFRA,QLAT,QLON)
        CALL FPSGCR (PLAT,PLON,QLAT,QLON,YFRA,ELAT,ELON)
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE FPSGCR (ALAT,ALON,BLAT,BLON,FRCT,FLAT,FLON)
C
C (FPSGCR = Find Point on Shortest Great Circle Route)
C
        REAL    ALAT,ALON,BLAT,BLON,FRCT,FLAT,FLON
C
C This routine, given the latitudes and longitudes of two points A and
C B on the surface of the globe and a fraction FRCT, between 0. and 1.,
C interpolates a point F along the shortest great circle route joining
C A to B such that the distance from A to F, divided by the distance
C from A to B, is equal to FRCT, and returns its latitude and longitude
C in FLAT and FLON.
C
C Define the constants used to convert from degrees to radians and
C vice-versa.
C
        DATA DTOR / .017453292519943 /
        DATA RTOD / 57.2957795130823 /
C
C Compute the X, Y, and Z coordinates (on a unit sphere) of the point B.
C
        CALL LLPXYZ (BLAT,BLON,XCPB,YCPB,ZCPB)
C
C Rotating about the Z axis by the angle -ALON would carry A into the XZ
C plane.
C
        CALL ROTXYZ (3,-ALON,XCPB,YCPB,ZCPB)
C
C Then, rotating about the Y axis by the angle ALAT would carry A into
C the point on the X axis with coordinates (1,0,0).
C
        CALL ROTXYZ (2,ALAT,XCPB,YCPB,ZCPB)
C
C Then, rotating about the X axis by the angle ALPH = -ATAN(ZCPB/YCPB)
C would leave the position of A unchanged but carry B into a point in
C the XY plane.
C
        IF (ZCPB.NE.0..OR.YCPB.NE.0.) THEN
          ALPH=-RTOD*ATAN2(ZCPB,YCPB)
        ELSE
          ALPH=0.
        END IF
C
        CALL ROTXYZ (1,ALPH,XCPB,YCPB,ZCPB)
C
C The angle BETA from A to B can now be computed easily.
C
        IF (XCPB.NE.0..OR.YCPB.NE.0.) THEN
          BETA=ATAN2(YCPB,XCPB)
        ELSE
          BETA=0.
        END IF
C
C Interpolate a point at the desired position between the points A and
C B, map it back to its original position on the great circle route from
C A to B, and get its latitude and longitude to return to the caller.
C
        GAMA=FRCT*BETA
C
        XCPF=COS(GAMA)
        YCPF=SIN(GAMA)
        ZCPF=0.
C
        CALL ROTXYZ (1,-ALPH,XCPF,YCPF,ZCPF)
        CALL ROTXYZ (2,-ALAT,XCPF,YCPF,ZCPF)
        CALL ROTXYZ (3, ALON,XCPF,YCPF,ZCPF)
C
        CALL XYZLLP (XCPF,YCPF,ZCPF,FLAT,FLON)
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE LLPXYZ (RLAT,RLON,RVOX,RVOY,RVOZ)
C
C (LLPXYZ = Lat/Lon Position to XYZ coordinates)
C
C Given the latitude and longitude of a point on the globe, return its
C X, Y, and Z coordinates.
C
        DATA DTOR / .017453292519943 /
C
        RVOX=COS(DTOR*RLAT)*COS(DTOR*RLON)
        RVOY=COS(DTOR*RLAT)*SIN(DTOR*RLON)
        RVOZ=SIN(DTOR*RLAT)
C
        RETURN
C
      END


      SUBROUTINE XYZLLP (RVOX,RVOY,RVOZ,RLAT,RLON)
C
C (XYZLLP = XYZ coordinates to Lat/Lon Position)
C
C Given the X, Y, and Z coordinates of a point on the globe, return its
C latitude and longitude.
C
        DATA RTOD / 57.2957795130823 /
C
        RLAT=RTOD*ASIN(RVOZ)
C
        IF (RVOX.NE.0.OR.RVOY.NE.0.) THEN
          RLON=RTOD*ATAN2(RVOY,RVOX)
        ELSE
          RLON=0.
        END IF
C
        RETURN
C
      END


      SUBROUTINE ROTXYZ (IAXS,ANGL,XCRD,YCRD,ZCRD)
C
C (ROTXYZ = ROTate a point defined by its X, Y, and Z coordinates)
C
C This is a modified version of a routine in the NCAR Graphics package,
C which is used in some of the examples.  It rotates the point having
C coordinates (XCRD,YCRD,ZCRD) by an angle ANGL about the axis specified
C by IAXS (1 for the X axis, 2 for the Y axis, 3 for the Z axis).
C
C One assumes a right-handed system with X, Y, and Z axes.  Rotating by
C an angle "a" about the X axis maps the point (x,y,z) into the point
C (x',y',z'), where
C
C       x' = x
C       y' = y cos(a) - z sin(a)
C       z' = z cos(a) + y sin(a)
C
C A positive value of "a" represents rotation in the direction from the
C Y axis to the Z axis.
C
C Similarly, rotating by an angle "a" about the Y axis maps the point
C (x,y,z) into the point (x',y',z'), where
C
C       x' = x cos(a) + z sin(a)
C       y' = y
C       z' = z cos(a) - x sin(a)
C
C A positive value of "a" represents rotation in the direction from the
C Z axis to the X axis.
C
C Rotating by an angle "a" about the Z axis maps the point (x,y,z) into
C the point (x',y',z'), where
C
C       x' = x cos(a) - y sin(a)
C       y' = y cos(a) + x sin(a)
C       y' = y
C
C A positive value of "a" represents rotation in the direction from the
C X axis to the Y axis.
C
C Define a multiplicative constant to convert from degrees to radians.
C
        DATA DTOR / .017453292519943 /
C
C Trigonometry.
C
        SINA=SIN(DTOR*ANGL)
        COSA=COS(DTOR*ANGL)
C
        XTMP=XCRD
        YTMP=YCRD
        ZTMP=ZCRD
C
        IF (IAXS.EQ.1) THEN
          YCRD=YTMP*COSA-ZTMP*SINA
          ZCRD=ZTMP*COSA+YTMP*SINA
        ELSE IF (IAXS.EQ.2) THEN
          XCRD=XTMP*COSA+ZTMP*SINA
          ZCRD=ZTMP*COSA-XTMP*SINA
        ELSE
          XCRD=XTMP*COSA-YTMP*SINA
          YCRD=YTMP*COSA+XTMP*SINA
        END IF
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
