*=======================================================================
*
* PGSBOX 4.8 - draw curvilinear coordinate axes for PGPLOT.
* Copyright (C) 1997-2011, Mark Calabretta
*
* This file is part of PGSBOX.
*
* PGSBOX is free software: you can redistribute it and/or modify it
* under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* PGSBOX is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with PGSBOX.  If not, see http://www.gnu.org/licenses.
*
* Correspondence concerning PGSBOX may be directed to:
*   Internet email: mcalabre@atnf.csiro.au
*   Postal address: Dr. Mark Calabretta
*                   Australia Telescope National Facility, CSIRO
*                   PO Box 76
*                   Epping NSW 1710
*                   AUSTRALIA
*
* Author: Mark Calabretta, Australia Telescope National Facility
* http://www.atnf.csiro.au/~mcalabre/index.html
* $Id: pgsbox.f,v 4.8.1.1 2011/08/15 08:07:07 cal103 Exp cal103 $
*=======================================================================
*
* PGSBOX draws and labels a curvilinear coordinate grid.  The caller
* must provide a separate external function, NLFUNC, to define the
* non-linear coordinate transformation.
*
* PGLBOX, a simplified ENTRY point to PGSBOX, has been provided for
* drawing simple linear axes without the need to specify NLFUNC.
* PGLBOX allows simplified access to formatting control for labelling
* world coordinate axes that is not provided by the standard PGPLOT
* routines, PGBOX or PGTBOX.  PGLBOX uses the world coordinate range
* set by a prior call to PGSWIN and omits the following arguments:
*
*   BLC, TRC, NLFUNC, NLC, NLI, NLD, NLCPRM, NLIPRM, NLDPRM
*
* The remaining arguments are specified in the same order as PGSBOX.
*
* Given:
*   BLC       R(2)      Cartesian coordinates of the bottom left-hand
*   TRC       R(2)      corner and top right-hand corner of the frame.
*                       Any convenient Cartesian system may be used in
*                       conjunction with the non-linear transformation
*                       function, NLFUNC, described below.
*
*                       For example, FITS images have pixel coordinate
*                       (1,1) at the centre of the pixel in the bottom
*                       left-hand corner.  Thus it would be
*                       appropriate to set BLC to (0.5,0.5) and TRC to
*                       (NAXIS1+0.5, NAXIS2+0.5).
*
*   IDENTS    C(3)*(*)  Identification strings:
*                         1: Name of the first world coordinate
*                            element used for axis labelling.
*                         2: Name of the second world coordinate.
*                         3: Title, written at top.
*
*   OPT       C(2)*(*)  Formatting control for the world coordinates,
*                       used for axis labelling (see notes 1 and 2):
*                         ' ': plain numeric
*                         'A': angle in degrees expressed in decimal
*                              notation normalized in the range [0,360).
*                         'B': as 'A' but normalized in the range
*                              (-180,180].
*                         'C': as 'A' but unnormalized.
*                         'D': angle in degrees expressed in sexagesimal
*                              notation, DD^MM'SS".S, normalized in the
*                              range [0,360).
*                         'E': as 'D' but normalized in the range
*                              (-180,180].
*                         'F': as 'D' but unnormalized.
*                         'G': angle in degrees expressed in hms
*                              notation, HHhMMmSSs.S, normalized in the
*                              range [0,24) hours.
*                         'H': as 'G' but normalized in the range
*                              (-12,12] hours.
*                         'I': as 'G' but unnormalized.
*                         'L': logarithmic (see note 2)
*                         'T': time in hours expressed as HH:MM:SS.S
*                         'Y': Modified Julian Date to be expressed as
*                              a Gregorian calendar date, YYYY/MM/DD.
*                              (MJD = JD - 2400000.5.)
*
*                       For the angular types, NLFUNC is assumed to
*                       return the angle in degrees whereupon it will
*                       be formatted for display in the specified way.
*                       For example, an angle of -417.2958 degrees
*                       (returned by NLFUNC):
*
*                         'A':  302^.7042
*                         'B':  -57^.2958
*                         'C': -417^.2958
*                         'D':  302^42'15"
*                         'E':  -57^17'45"
*                         'F': -417^17'45"
*                         'H':   20h10m49s
*                         'I':   -3h49m11s
*                         'J':  -27h49m11s
*
*                       These are properly superscripted.
*
*   LABCTL    I         Decimal encoded grid labelling control:
*                            -1: Accumulate information on grid labels
*                                (see note 3) but do not write them.
*                             0: Let PGSBOX decide what edges to label.
*                             1: Label bottom of frame with the first
*                                world coordinate.
*                             2: Label bottom of frame with the second
*                                world coordinate.
*                            10: ... left side of frame.
*                            20: ... left side of frame.
*                           100: ... top of frame.
*                           200: ... top of frame.
*                          1000: ... right side of frame.
*                          2000: ... right side of frame.
*                         10000: Write labels from information
*                                accumulated in previous calls without
*                                drawing grid lines or tick marks.
*
*                       LABCTL = 0 usually gets what you want.
*                       LABCTL = 3333 labels all sides with both world
*                       coordinates.
*
*   LABDEN    I         Decimal encoded labelling density control for
*                       use where PGSBOX is called upon to determine a
*                       suitable grid spacing (e.g. via NG1 = 0,
*                       GRID1(0) = 0).  LABDEN = 100*D2 + D1 where
*                       D1, and D2 are the approximate number of grid
*                       lines for the first and second world
*                       coordinate.  LABDEN = 0 is effectively the same
*                       as LABDEN = 808.
*
*   CI        I(7)      Table of predefined colours established by
*                       calls to PGSCR.  This is used to control the
*                       colour used for different parts of the plot.
*                       CI table entries are used as follows:
*
*                                                world
*                                              coordinate
*                       Index      usage        element
*                       -----  --------------  ----------
*                         1      grid lines        1
*                         2      grid lines        2
*                         3    numeric labels      1
*                         4    numeric labels      2
*                         5    axis annotation     1
*                         6    axis annotation     2
*                         7        title           -
*
*                       For example, CI(3) is used for numeric labels
*                       for the first world coordinate.
*
*                       Colour selection is disabled for component J
*                       if CI(J) < 0.
*
*   GCODE(2)  I         Code for the type of grid to draw for each
*                       world coordinate:
*                         0: No grid or tick marks.
*                         1: Tick marks (on all edges).
*                         2: Full coordinate grid.
*
*                       Tick marks can be restricted to particular
*                       edges of the frame; a negative GCODE is
*                       interpreted as a decimal encoded control
*                       variable:
*                            -1: bottom
*                           -10: left
*                          -100: top
*                         -1000: right
*
*                       The digit scales the basic tick length.  For
*                       example, GCODE(1) = -102 restricts tick marks
*                       for the first world coordinate to the bottom
*                       and top edges of the frame.  Those on the
*                       bottom will be twice the length specified by
*                       TIKLEN.
*
*   TIKLEN    D         Tick length, in mm.  Negative values produce
*                       outside tick marks.
*
*   NG1       I         Upper array index for GRID1.
*
*   GRID1     D(0:NG1)  Grid values in WORLD(1) in the same units as
*                       returned by NLFUNC.
*
*                       If NG1 is zero, then
*                         a: if GRID1(0) is greater than zero it
*                            defines a uniform grid spacing.
*                         b: if GRID1(0) is zero a suitable spacing
*                            will be determined (see LABDEN).
*                         c: if GRID1(0) is less than zero then no
*                            grid lines will be drawn.
*
*                       If NG1 is greater than zero, then GRID1(0) is
*                       ignored.
*
*   NG2       I         Upper array index for GRID2.
*
*   GRID2     D(0:NG2)  Grid values in WORLD(2) in the same units as
*                       returned by NLFUNC, interpreted the same way
*                       as GRID1.
*
*   DOEQ      L         If NG1 = NG2 = 0, and GRID1(0) = 0D0 and/or
*                       GRID2(0) = 0D0, then choose the same grid
*                       spacing for each world coordinate.
*
*   NLFUNC    Ext       Non-linear coordinate function, see below.
*
*   NLC       I         Number of elements in NLCPRM (must be >0).
*
*   NLI       I         Number of elements in NLIPRM (must be >0).
*
*   NLD       I         Number of elements in NLDPRM (must be >0).
*
* Given and/or returned:
*   NLCPRM    C(NLC)*1  Character coefficients for NLFUNC.
*
*   NLIPRM    I(NLI)    Integer coefficients for NLFUNC.
*
*   NLDPRM    D(NLD)    Double precision coefficients for NLFUNC.
*
*   NC        I         Upper array index for CACHE (see note 3).
*
*   IC        I         Current number of entries in the CACHE table.
*                       Should be set to -1 on the first call to
*                       PGSBOX (see note 3).
*
*   CACHE     D(4,0:NC) Table of points where the tick marks or grid
*                       lines cross the frame (see note 3).
*                         1: Frame segment
*                            1: bottom
*                            2: left
*                            3: top
*                            4: right
*                         2: X or Y-Cartesian coordinate.
*                         3: World coordinate element (1 or 2).
*                         4: Value.
*
*                       CACHE(,0) is used to cache the extrema of the
*                       coordinate values between calls.  CACHE(1,NC-1)
*                       is also used to store related information.
*
*                       CACHE(,NC) will contain the margin widths in
*                       Cartesian coordinates when the labels are
*                       produced (i.e. the same Cartesian system used
*                       for BLC and TRC).
*
* Returned:
*   IERR      I         Status return value:
*                         0: Success
*                         1: Initialization error
*                         2: Invalid coordinate system
*                         3: Cache overflow (see note 3).
*
* Notes:
*   1: Where a logarithmic world coordinate type is indicated PGSBOX
*      chooses grid lines and labels on the basis that the value
*      returned by NLFUNC is a base 10 logarithm.  PGSBOX does not
*      itself take logarithms or antilogarithms.  For example, if the
*      range of values returned by NLFUNC were 0.9 - 2.5, then PGSBOX
*      would draw a subset of the following set of grid lines and
*      labels:
*
*                 value          label
*          ------------------    -----
*          0.9031 = log10(8)      8
*          0.9542 = log10(9)      9
*          1.0000 = log10(10)    10**1
*          1.3010 = log10(20)     2
*          1.4771 = log10(30)     3
*          1.6021 = log10(40)     4
*          1.6990 = log10(50)     5
*          1.7782 = log10(60)     6
*          1.8451 = log10(70)     7
*          1.9031 = log10(80)     8
*          1.9542 = log10(90)     9
*          2.0000 = log10(100)   10**2
*          2.3010 = log10(200)    2
*          2.4771 = log10(300)    3
*
*      The subset chosen depends on the coordinate increment as
*      specified by the caller (e.g. via NG1 = 0, GRID1(0) > 0) or as
*      deduced by PGSBOX from the required density of grid lines
*      specified in the LABDEN argument.  The selection is made
*      according to the following table:
*
*          increment          grid lines
*         (0.00,0.12]   1, 2, 3, 4, 5, 6, 7, 8, 9
*         (0.12,0.18]   1, 2, 3, 4, 5,    7
*         (0.18,0.23]   1, 2, 3,    5,    7
*         (0.23,0.28]   1, 2, 3,    5
*         (0.28,0.40]   1, 2,       5
*         (0.40,0.70]   1, 3
*         (0.70,1.00]   1
*
*      For increments greater than 1 the nearest integer is used.
*
*   2: PGSBOX will attempt to handle discontinuities in angle, such as
*      may occur when cycling through 360 degrees, wherever the
*      discontinuity may occur, for example
*
*        359 -> 0 -> 1
*
*      or
*
*        179 -> 180 -> -179
*
*      Only single cycles are detected, so the sequence
*
*        -360 -> 0 -> ... -> 359 -> 0 -> ... -> 359
*
*      would not be handled properly.  In such cases NLFUNC should be
*      changed to return a normalized angle, or else a continuous
*      sequence.
*
*   3: PGSBOX maintains a table of axis crossings, CACHE, in which it
*      stores information used for axis labelling.  The caller need not
*      normally be concerned about the use of this table other than to
*      provide sufficient space.  Typically, NC = 256 should be enough;
*      if not, IERR = 3 will be returned.
*
*      However, a coordinate grid may be produced via multiple calls to
*      PGSBOX with deferment of axis labelling.  This might be done to
*      change the pen colour and/or thickness for different sets of
*      grid lines.  The table accumulates information from each
*      successive call until the labels are produced.  The table index,
*      IC, may be reset to zero by the caller to discard the
*      information collected.
*
*      The extrema of the world coordinate elements are stored in
*      CACHE(,0).  When a coordinate grid is plotted with multiple
*      calls to PGSBOX the initial call should always have IC set to -1
*      to signal that PGSBOX needs to determine the extrema.  On
*      subsequent calls with IC non-negative PGSBOX uses the extrema
*      cached from the first call.  This can speed up execution
*      considerably.
*
*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*
* The curvilinear coordinate grid is defined by external function
* NLFUNC whose interface is defined as follows:
*
* Given:
*   OPCODE    I         Transformation code:
*                         +2: Compute a set of Cartesian coordinates
*                             that describe a path between this and
*                             the previous pair of world coordinates
*                             remembered from the last call with
*                             OPCODE = +1 or +2.  Usually only takes
*                             a single step unless traversing a
*                             discontinuity or some other
*                             irregularity (see explanation below).
*                         +1: Compute Cartesian coordinates from world
*                             coordinates.
*                          0: Initialize.
*                         -1: Compute world coordinates from Cartesian
*                             coordinates.
*
*                       N.B. NLFUNC must not change the input
*                       coordinates; that is the world coordinates for
*                       OPCODEs = +2 and +1, or the Cartesian
*                       coordinates for OPCODE = -1.
*
*   NLC       I         Number of elements in NLCPRM (must be >0).
*
*   NLI       I         Number of elements in NLIPRM (must be >0).
*
*   NLD       I         Number of elements in NLDPRM (must be >0).
*
* Given and/or returned:
*   NLCPRM    C(NLC)*1  Character array.
*
*   NLIPRM    I(NLI)    Integer coefficients.
*
*   NLDPRM    D(NLD)    Double precision coefficients.
*
*   WORLD     D(2)      World coordinates.  Given if OPCODE > 0,
*                       returned if OPCODE < 0.
*
*   XY        D(2)      Cartesian coordinates.  Given if OPCODE < 0,
*                       returned if OPCODE > 0.
*
*   CONTRL    I         Control flag for OPCODE = +2:
*                         0: Normal state.
*                         1: Force PGSBOX to flush its plotting buffer
*                            and call NLFUNC again with the same world
*                            coordinates.
*                         2: Force PGSBOX to call NLFUNC again with
*                            the same world coordinates (without
*                            flushing its plotting buffer).
*
*   CONTXT    D(20)     Context elements for OPCODE = +2.
*
* Returned:
*   IERR      I         Status return value:
*                         0: Success.
*                         1: Invalid parameters.
*                         2: Invalid world coordinate.
*                         3: Invalid Cartesian coordinate.
*
*                       The following status returns are recognized for
*                       opcodes +2 and +1
*                         -1: Accept the returned (x,y) coordinates but
*                             do not consider this as one end of a
*                             crossing segment for labelling world
*                             coordinate 1.
*                         -2: Ditto for world coordinate 2.
*                         -3: Ditto for world coordinates 1 and 2.
*
* PGSBOX passes its NLCPRM, NLIPRM, and NLDPRM adjustable size array
* arguments of length NLC, NLI, and NLD to NLFUNC without
* modification.  Comments within NLFUNC should specify the parameters
* it wants passed to it via these arrays.
*
* PGSBOX first calls NLFUNC with OPCODE = 0 to cause it to initialize
* its work arrays (if necessary).  It then uses OPCODE = -1 to
* determine the range of world coordinate values.  It anchors the
* start of each coordinate grid line with a call with OPCODE = +1, and
* then tracks it with OPCODE = +2.
*
* The CONTXT array is also passed to NLFUNC without modification to
* allow it to preserve state information between calls for OPCODE = 2.
* In particular, NLFUNC can use this to detect discontinuities in the
* grid lines.
*
* The CONTRL argument is provided so that NLFUNC can force PGSBOX to
* call it again with or without flushing its plotting buffer.  This
* may be needed when plotting a grid line through a discontinuity.
* PGSBOX does not modify CONTRL.
*
* Notes:
*   1: NLFUNC must not change the input coordinates; that is the world
*      coordinates for OPCODEs = +1 and +2, or the Cartesian coordinates
*      for OPCODE = -1.
*
*   2: NLFUNC must define a single-valued function, that is, each
*      Cartesian coordinate (x,y) must map to a unique world coordinate
*      pair (xi,eta).
*
*   3: Notwithstanding the fact that PGSBOX declares NLCPRM, NLIPRM,
*      and NLDPRM as single dimension arrays of length NLC, NLI, and
*      NLD, NLFUNC may treat these as higher-dimensional arrays, for
*      example, NLDPRM(2,NLD).  (The FORTRAN standard requires that
*      only the last dimension is adjustable.)
*
*=======================================================================
      SUBROUTINE PGSBOX (BLC_, TRC_, IDENTS, OPT, LABCTL, LABDEN, CI,
     :  GCODE, TIKLEN, NG1, GRID1, NG2, GRID2, DOEQ, NLFUNC, NLC, NLI,
     :  NLD, NLCPRM, NLIPRM, NLDPRM, NC, IC, CACHE, IERR)
*-----------------------------------------------------------------------
      INTEGER   BUFSIZ
      PARAMETER (BUFSIZ = 2048)

      LOGICAL   DOEDGE, DOEQ, DOLBOX, FULLSM, GETEND, INSIDE, ISANGL(2),
     :          LABLOK, MAJOR, OVERFL, PREVIN
      INTEGER   CI(7), CI0, CJ(7), CONTRL, DENS(2), FSEG, GCODE(2), IC,
     :          ID, IERR, IM, ISTEP, IW0, IWJ, IWK, IX, IY, IYPREV, J,
     :          K, KX, L, LABCTL, LABDEN, LDIV(2), LTABL(6,2:6), NC,
     :          NG(2), NG1, NG2, NLC, NLD, NLI, NLIPRM(NLI), NP,
     :          NSTEP(2), NWJ, NX, NY, TCODE(2,4)
      REAL      BLC(2), BLC_(2), S, TRC(2), TRC_(2), WXY(4), X1, X2,
     :          XPOINT, XR(BUFSIZ), XSCL, XSPAN, XTOL, XVP1, XVP2,
     :          Y1, Y2, YR(BUFSIZ), YSCL, YSPAN, YTOL, YVP1, YVP2
      DOUBLE PRECISION CONTXT(20), CACHE(4,0:NC), DW(2), DX, DY, FACT,
     :          G0(2), GSTEP(2), GRID1(0:NG1), GRID2(0:NG2),
     :          NLDPRM(NLD), STEP, SW(2), TIKLEN, TMP, VMAX(2,2),
     :          VMIN(2,2), W1PREV, W1X0, W2PREV, W2X0, WJUMP, WMAX(2),
     :          WMIN(2), WORLD(9), XY(9)
      CHARACTER FTYPE(2), IDENTS(3)*(*), NLCPRM(NLC)*1, OPT(2)*(*)

      EXTERNAL NLFUNC

*     Approximate number of grid lines for each coordinate.
      INTEGER DENS0
      PARAMETER (DENS0 = 8)

*     Double precision round-off tolerance.
      DOUBLE PRECISION TOL
      PARAMETER (TOL = 1D-8)

*     Number of steps per grid line.
      DATA NSTEP /80, 80/

*     Table of logarithmic grid values.
      DATA LTABL /3, 10,  0,  0,  0,  0,
     :            2,  5, 10,  0,  0,  0,
     :            2,  3,  5, 10,  0,  0,
     :            2,  3,  5,  7, 10,  0,
     :            2,  3,  4,  5,  7, 10/

*     These are to stop compiler messages about uninitialized variables.
      DATA IW0 /0/
      DATA LABLOK, PREVIN /2 * .FALSE./
      DATA G0 /2 * 0D0/
      DATA W1X0, W1PREV, W2X0, W2PREV /4 * 0D0/
      DATA WORLD, XY /9 * 0D0, 9 * 0D0/
*-----------------------------------------------------------------------
*  Initialize.
      DOLBOX = .FALSE.

      CALL NLFUNC (0, NLC, NLI, NLD, NLCPRM, NLIPRM, NLDPRM, WORLD,
     :  XY, CONTRL, CONTXT, IERR)
*     Quick return for now.
      IF (IERR.NE.0) THEN
        IERR = 1
        RETURN
      END IF

      BLC(1) = BLC_(1)
      BLC(2) = BLC_(2)
      TRC(1) = TRC_(1)
      TRC(2) = TRC_(2)

      DOEDGE = GCODE(1).NE.2 .AND. GCODE(2).NE.2

*-----------------------------------------------------------------------
      GO TO 10

      ENTRY PGLBOX (IDENTS, OPT, LABCTL, LABDEN, CI, GCODE, TIKLEN, NG1,
     :  GRID1, NG2, GRID2, DOEQ, NC, IC, CACHE, IERR)

      DOLBOX = .TRUE.

      BLC(1) = 0.0
      BLC(2) = 0.0
      TRC(1) = 1.0
      TRC(2) = 1.0

      DOEDGE = .TRUE.

      CONTRL = 0
      IERR = 0

 10   CONTINUE

*-----------------------------------------------------------------------
      IF (NC.LT.1) THEN
        IERR = 3
        RETURN
      END IF

      NG(1) = NG1
      NG(2) = NG2

      FTYPE(1) = OPT(1)(1:1)
      FTYPE(2) = OPT(2)(1:1)

*     Extend the PGPLOT window and rescale it.
      CALL PGQVP (0, XVP1, XVP2, YVP1, YVP2)
      CALL PGQWIN (WXY(1), WXY(2), WXY(3), WXY(4))
      CALL PGSVP (0.0, 1.0, 0.0, 1.0)
      XSCL = (TRC(1)-BLC(1))/(XVP2-XVP1)
      YSCL = (TRC(2)-BLC(2))/(YVP2-YVP1)
      CALL PGSWIN (BLC(1)-XSCL*XVP1, TRC(1)+XSCL*(1.0-XVP2),
     :             BLC(2)-YSCL*YVP1, TRC(2)+YSCL*(1.0-YVP2))

*     Determine initial colour and set colour indices.
      CALL PGQCI (CI0)
      DO 20 J = 1, 7
        IF (CI(J).GE.0) THEN
          CJ(J) = CI(J)
        ELSE
          CJ(J) = CI0
        END IF
 20   CONTINUE

*     Labels only?
      IF (LABCTL.GE.10000) GO TO 130


      XSPAN = WXY(2) - WXY(1)
      YSPAN = WXY(4) - WXY(3)
      XTOL  = XSPAN*TOL
      YTOL  = YSPAN*TOL

*  Find world coordinate ranges:
*  WMIN(1:2)  ...lower limit of each world coordinate element.
*  WMAX(1:2)  ...upper limit of each world coordinate element.
*  GSTEP(1:2) ...spacing between grid lines for each element.
*  DW(1:2)    ...span from WMIN to WMAX.
*  SW(1:2)    ...increment in world coordinate between plotting steps.
*  NSTEP(1:2) ...number of plotting steps between WMIN and WMAX.
      FULLSM = .FALSE.
      IF (IC.GE.0 .AND. IC.LT.NC-1) FULLSM = CACHE(1,NC-1).EQ.1D0

      IF (IC.GE.0 .AND. (FULLSM .OR. DOEDGE)) THEN
*       Use extrema cached from a previous call.
        WMIN(1) = CACHE(1,0)
        WMAX(1) = CACHE(2,0)
        WMIN(2) = CACHE(3,0)
        WMAX(2) = CACHE(4,0)

      ELSE
*       Do a coarse search to find approximate ranges.
        WMIN(1) =  1D99
        WMAX(1) = -1D99
        WMIN(2) =  1D99
        WMAX(2) = -1D99

*       Need to consider cycles in angle through 360 degrees.
        ISANGL(1) = INDEX('ABCDEFGHI',FTYPE(1)).NE.0
        ISANGL(2) = INDEX('ABCDEFGHI',FTYPE(2)).NE.0
        VMIN(1,1) =  1D99
        VMIN(1,2) =  1D99
        VMAX(1,1) = -1D99
        VMAX(1,2) = -1D99
        VMIN(2,1) =  1D99
        VMIN(2,2) =  1D99
        VMAX(2,1) = -1D99
        VMAX(2,2) = -1D99

*       Sample coordinates on a 50 x 50 grid.
        NX = 50
        NY = 50
        DX = DBLE(TRC(1)-BLC(1))/NX
        DY = DBLE(TRC(2)-BLC(2))/NY

        K = 0
        IYPREV = -1
        DO 40 IY = 0, NY
          XY(2) = BLC(2) + IY*DY

*         Sample the edges only?
          KX = 1
          IF (DOEDGE) THEN
            IF (IY.NE.0 .AND. IY.NE.NY) KX = NX
          END IF

          DO 30 IX = 0, NX, KX
            XY(1) = BLC(1) + IX*DX

            IF (DOLBOX) THEN
              WORLD(1) = WXY(1) + XY(1)*XSPAN
              WORLD(2) = WXY(3) + XY(2)*YSPAN
            ELSE
              CALL NLFUNC (-1, NLC, NLI, NLD, NLCPRM, NLIPRM,
     :          NLDPRM, WORLD, XY, CONTRL, CONTXT, IERR)
            END IF

            IF (IERR.EQ.0) THEN
              K = K + 1

              IF (ISANGL(1)) THEN
                IF (K.EQ.1) W1X0 = WORLD(1)

                IF (IY.NE.IYPREV) THEN
                  W1PREV = W1X0
                  W1X0 = WORLD(1)
                END IF

*               Iron out jumps.
                WJUMP = WORLD(1) - W1PREV
                IF (ABS(WJUMP).GT.180D0) THEN
                  WJUMP = WJUMP + SIGN(180D0,WJUMP)
                  WJUMP = 360D0*INT(WJUMP/360D0)
                  WORLD(1) = WORLD(1) - WJUMP
                END IF

                W1PREV = WORLD(1)
                IYPREV = IY
              END IF

              IF (ISANGL(2)) THEN
                IF (K.EQ.1) W2X0 = WORLD(2)

                IF (IY.NE.IYPREV) THEN
                  W2PREV = W2X0
                  W2X0 = WORLD(2)
                END IF

*               Iron out jumps.
                WJUMP = WORLD(2) - W2PREV
                IF (ABS(WJUMP).GT.180D0) THEN
                  WJUMP = WJUMP + SIGN(180D0,WJUMP)
                  WJUMP = 360D0*INT(WJUMP/360D0)
                  WORLD(2) = WORLD(2) - WJUMP
                END IF

                W2PREV = WORLD(2)
                IYPREV = IY
              END IF

              IF (WORLD(1).LT.WMIN(1)) WMIN(1) = WORLD(1)
              IF (WORLD(1).GT.WMAX(1)) WMAX(1) = WORLD(1)
              IF (WORLD(2).LT.WMIN(2)) WMIN(2) = WORLD(2)
              IF (WORLD(2).GT.WMAX(2)) WMAX(2) = WORLD(2)

              IF (ISANGL(1)) THEN
*               Normalize to the range [0,360).
                WORLD(1) = MOD(WORLD(1), 360D0)
                IF (WORLD(1).LT.0D0) WORLD(1) = WORLD(1) + 360D0
                IF (WORLD(1).LT.VMIN(1,1)) VMIN(1,1) = WORLD(1)
                IF (WORLD(1).GT.VMAX(1,1)) VMAX(1,1) = WORLD(1)

*               Normalize to the range (-180,180].
                IF (WORLD(1).GT.180D0) WORLD(1) = WORLD(1) - 360D0
                IF (WORLD(1).LT.VMIN(1,2)) VMIN(1,2) = WORLD(1)
                IF (WORLD(1).GT.VMAX(1,2)) VMAX(1,2) = WORLD(1)
              END IF

              IF (ISANGL(2)) THEN
*               Normalize to the range [0,360).
                WORLD(2) = MOD(WORLD(2), 360D0)
                IF (WORLD(2).LT.0D0) WORLD(2) = WORLD(2) + 360D0
                IF (WORLD(2).LT.VMIN(2,1)) VMIN(2,1) = WORLD(2)
                IF (WORLD(2).GT.VMAX(2,1)) VMAX(2,1) = WORLD(2)

*               Normalize to the range (-180,180].
                IF (WORLD(2).GT.180D0) WORLD(2) = WORLD(2) - 360D0
                IF (WORLD(2).LT.VMIN(2,2)) VMIN(2,2) = WORLD(2)
                IF (WORLD(2).GT.VMAX(2,2)) VMAX(2,2) = WORLD(2)
              END IF
            END IF
 30       CONTINUE
 40     CONTINUE

        IF (K.EQ.0) THEN
*         No valid coordinates found within the frame.
          IERR = 2
          GO TO 999
        END IF

*       Check for cycles in angle.
        DO 50 J = 1, 2
          IF (ISANGL(J)) THEN
            IF (WMAX(J)-WMIN(J).LT.360D0 .AND.
     :          WMAX(J)-WMIN(J).GT.VMAX(J,1)-VMIN(J,1)+TOL) THEN
*             Must have a cycle, preserve the sign.
              IF (WMAX(J).GE.0D0) THEN
                WMIN(J) = VMIN(J,1)
                WMAX(J) = VMAX(J,1)
              ELSE
                WMIN(J) = VMIN(J,1) - 360D0
                WMAX(J) = VMAX(J,1) - 360D0
              END IF
            END IF

            IF (WMAX(J)-WMIN(J).LT.360D0 .AND.
     :          WMAX(J)-WMIN(J).GT.VMAX(J,2)-VMIN(J,2)+TOL) THEN
*             Must have a cycle, preserve the sign.
              IF (WMAX(J).GE.0D0) THEN
                IF (VMAX(J,2).GE.0D0) THEN
                  WMIN(J) = VMIN(J,2)
                  WMAX(J) = VMAX(J,2)
                ELSE
                  WMIN(J) = VMIN(J,2) + 360D0
                  WMAX(J) = VMAX(J,2) + 360D0
                END IF
              ELSE
                IF (VMAX(J,2).LT.0D0) THEN
                  WMIN(J) = VMIN(J,2)
                  WMAX(J) = VMAX(J,2)
                ELSE
                  WMIN(J) = VMIN(J,2) - 360D0
                  WMAX(J) = VMAX(J,2) - 360D0
                END IF
              END IF
            END IF
          END IF
 50     CONTINUE

*       Cache extrema for subsequent calls.
        CACHE(1,0) = WMIN(1)
        CACHE(2,0) = WMAX(1)
        CACHE(3,0) = WMIN(2)
        CACHE(4,0) = WMAX(2)

*       Was full sampling done?
        IF (DOEDGE) THEN
          CACHE(1,NC-1) = 0D0
        ELSE
          CACHE(1,NC-1) = 1D0
        END IF

        IC = 0
      END IF

*     Choose an appropriate grid spacing.
      IF (LABDEN.GT.0) THEN
*       User specified grid density.
        DENS(1) = MOD(LABDEN,100)
        DENS(2) = LABDEN/100
        IF (DENS(1).EQ.0) DENS(1) = DENS0
        IF (DENS(2).EQ.0) DENS(2) = DENS0
      ELSE
*       Default grid density.
        DENS(1) = DENS0
        DENS(2) = DENS0
      END IF

      IF (NG(1).EQ.0) G0(1) = GRID1(0)
      IF (NG(2).EQ.0) G0(2) = GRID2(0)

      DO 60 J = 1, 2
        IF (J.EQ.1) THEN
          K = 2
        ELSE
          K = 1
        END IF

        IF (NG(J).EQ.0 .AND. G0(J).LT.0D0) THEN
*         Defeat grid lines and tick marks.
          GCODE(J) = 0
        END IF

        IF (NG(J).EQ.0 .AND. G0(J).GT.0D0) THEN
*         User-specified, directly.
          GSTEP(J) = G0(J)
        ELSE IF (DOEQ .AND. NG(K).EQ.0 .AND. G0(K).GT.0D0) THEN
*         User-specified, indirectly.
          GSTEP(J) = G0(K)
        ELSE
*         Left to us to choose.  Even if grid lines are not drawn for
*         coordinate J, i.e. NG(J) = 0 and G0(J) < 0, GSTEP(J) is still
*         needed because it is used to deduce SW(J) which is used for
*         drawing grid lines for the other coordinate.
          DW(J) = WMAX(J) - WMIN(J)
          STEP = DW(J)/DENS(J)

          FACT = 1D0
          IF (INDEX('GHI',FTYPE(J)).NE.0) THEN
*           Rescale degrees to hours.
            FACT = 1D0/15D0
            STEP = STEP*FACT
          ELSE IF (FTYPE(J).EQ.'Y' .AND. STEP.LT.0.5D0) THEN
*           Calendar increment of less than 12h; use time format.
            FTYPE(J) = 'y'

*           Rescale days to hours.
            FACT = 24D0
            STEP = STEP*FACT
          END IF

          IF (INDEX('ABCDEF',FTYPE(J)).NE.0 .AND. STEP.GE.1D0) THEN
*           Angle with multi-degree increment.
            IF (STEP.LT.1.5D0) THEN
              STEP = 1D0
            ELSE IF (STEP.LT.3D0) THEN
              STEP = 2D0
            ELSE IF (STEP.LT.7D0) THEN
              STEP = 5D0
            ELSE IF (STEP.LT.12D0) THEN
              STEP = 10D0
            ELSE IF (STEP.LT.20D0) THEN
              STEP = 15D0
            ELSE IF (STEP.LT.40D0) THEN
              STEP = 30D0
            ELSE IF (STEP.LT.70D0) THEN
              STEP = 45D0
            ELSE IF (STEP.LT.120D0) THEN
              STEP = 90D0
            ELSE IF (STEP.LT.270D0) THEN
              STEP = 180D0
            ELSE IF (STEP.LT.520D0) THEN
              STEP = 360D0
            ELSE
              STEP = 360D0*INT(STEP/360D0 + 0.5)
            END IF

          ELSE IF (INDEX('GHITy',FTYPE(J)).NE.0 .AND.
     :      STEP.GE.1D0) THEN
*           Angle or time in hms format with multi-hour increment.
            IF (STEP.LT.1.5D0) THEN
              STEP = 1D0
            ELSE IF (STEP.LT.2.5D0) THEN
              STEP = 2D0
            ELSE IF (STEP.LT.3.5D0) THEN
              STEP = 3D0
            ELSE IF (STEP.LT.5D0) THEN
              STEP = 4D0
            ELSE IF (STEP.LT.7D0) THEN
              STEP = 6D0
            ELSE IF (STEP.LT.10D0) THEN
              STEP = 8D0
            ELSE IF (STEP.LT.15D0) THEN
              STEP = 12D0
            ELSE IF (STEP.LT.21D0) THEN
              STEP = 18D0
            ELSE IF (STEP.LT.36D0) THEN
              STEP = 24D0
            ELSE
              STEP = 24D0*INT(STEP/24D0 + 0.5)
            END IF

            STEP = STEP/FACT

          ELSE IF (INDEX('DEFGHITy',FTYPE(J)).NE.0 .AND.
     :      STEP.LT.1D0) THEN
*           Angle or time in sexagesimal format with sub-degree/hour
*           increment.
            FACT = FACT*60D0
            STEP = STEP*60D0
            IF (STEP.LT.1D0) THEN
*             Sub-minute increment.
              FACT = FACT*60D0
              STEP = STEP*60D0
            END IF

            IF (STEP.LT.1D0) THEN
*             Sub-second increment.
              TMP = 10D0**(INT(LOG10(STEP))-1)
              IF (1.5*TMP.GE.STEP) THEN
                STEP = TMP
              ELSE
                IF (3D0*TMP.GE.STEP) THEN
                  STEP = 2D0*TMP
                ELSE
                  IF (7D0*TMP.GE.STEP) THEN
                    STEP = 5D0*TMP
                  ELSE
                    STEP = 10D0*TMP
                  END IF
                END IF
              END IF

            ELSE
              IF (STEP.LT.1.5D0) THEN
                STEP = 1D0
              ELSE IF (STEP.LT.2.5D0) THEN
                STEP = 2D0
              ELSE IF (STEP.LT.3.5D0) THEN
                STEP = 3D0
              ELSE IF (STEP.LT.4.5D0) THEN
                STEP = 4D0
              ELSE IF (STEP.LT.5.5D0) THEN
                STEP = 5D0
              ELSE IF (STEP.LT.8D0) THEN
                STEP = 6D0
              ELSE IF (STEP.LT.11D0) THEN
                STEP = 10D0
              ELSE IF (STEP.LT.14D0) THEN
                STEP = 12D0
              ELSE IF (STEP.LT.18D0) THEN
                STEP = 15D0
              ELSE IF (STEP.LT.25D0) THEN
                STEP = 20D0
              ELSE IF (STEP.LT.45D0) THEN
                STEP = 30D0
              ELSE
                STEP = 60D0
              END IF
            END IF

            STEP = STEP/FACT

          ELSE IF (FTYPE(J).EQ.'Y') THEN
*           Calendar axis: use coded steps.
            IF (STEP.LT.15D0) THEN
*             Timespan of a few months; use multi-day increments.
              STEP = ANINT(STEP)
              IF (STEP.LT.1D0) THEN
                STEP = 1D0
              ELSE IF (STEP.GT.9D0) THEN
*               Fortnightly.
                STEP = 14D0
              ELSE IF (STEP.GT.4D0) THEN
*               Weekly.
                STEP = 7D0
              END IF

            ELSE IF (STEP.LT.270D0) THEN
*             Timespan of a few years; use multi-month increments.
              STEP = ANINT(STEP/30.44D0)
              IF (STEP.LT.1.5D0) THEN
                STEP = 1D0
              ELSE IF (STEP.LT.2.5D0) THEN
                STEP = 2D0
              ELSE IF (STEP.LT.3.5D0) THEN
                STEP = 3D0
              ELSE IF (STEP.LT.4.5D0) THEN
                STEP = 4D0
              ELSE
                STEP = 6D0
              END IF

*             Coding for multi-month increments.
              STEP = 100D0*STEP

            ELSE
*             Multi-year increments.
              STEP = ANINT(DW(J)/DENS(J)/365.25D0)
              IF (STEP.LT.1D0) THEN
                STEP = 1D0
              ELSE
                TMP = 10D0**INT(LOG10(STEP))

                IF (1.5D0*TMP.GE.STEP) THEN
                  STEP = TMP
                ELSE
                  IF (3D0*TMP.GE.STEP) THEN
                    STEP = 2D0*TMP
                  ELSE
                    IF (7D0*TMP.GE.STEP) THEN
                      STEP = 5D0*TMP
                    ELSE
                      STEP = 10D0*TMP
                    END IF
                  END IF
                END IF
              END IF

*             Coding for multi-year increments.
              STEP = 10000D0*STEP
            END IF

          ELSE
*           Just numbers.
            TMP = 10D0**INT(LOG10(STEP))
            IF (STEP.LT.1D0) TMP = TMP/10D0

            IF (1.5D0*TMP.GE.STEP) THEN
              STEP = TMP
            ELSE
              IF (3D0*TMP.GE.STEP) THEN
                STEP = 2D0*TMP
              ELSE
                IF (7D0*TMP.GE.STEP) THEN
                  STEP = 5D0*TMP
                ELSE
                  STEP = 10D0*TMP
                END IF
              END IF
            END IF

*           Adjust the step size for logarithmic values.
            IF (FTYPE(J).EQ.'L') THEN
              IF (STEP.GT.0.7D0) THEN
                LDIV(J) = 1
                STEP = NINT(STEP)
              ELSE
                IF (STEP.GT.0.4D0) THEN
                  LDIV(J) = 2
                ELSE IF (STEP.GT.0.28D0) THEN
                  LDIV(J) = 3
                ELSE IF (STEP.GT.0.23D0) THEN
                  LDIV(J) = 4
                ELSE IF (STEP.GT.0.18D0) THEN
                  LDIV(J) = 5
                ELSE IF (STEP.GT.0.12D0) THEN
                  LDIV(J) = 6
                ELSE
                  LDIV(J) = 9
                END IF

                STEP = 1D0/LDIV(J)
              END IF
            END IF
          END IF

          GSTEP(J) = STEP
        END IF
 60   CONTINUE

*     Equal grid spacing?
      IF (DOEQ .AND. NG(1).EQ.0 .AND. NG(2).EQ.0) THEN
        IF (GRID1(0).EQ.0D0 .AND. GRID2(0).EQ.0D0) THEN
          GSTEP(1) = MIN(GSTEP(1), GSTEP(2))
          GSTEP(2) = GSTEP(1)
        ELSE IF (GRID1(0).EQ.0D0) THEN
          GSTEP(1) = GSTEP(2)
        ELSE IF (GRID2(0).EQ.0D0) THEN
          GSTEP(2) = GSTEP(1)
        END IF
      END IF

*     Fine tune the end points.
      DO 70 J = 1, 2
        IF (FTYPE(J).EQ.'L') THEN
          WMIN(J) = AINT(WMIN(J)-1D0)
          WMAX(J) = AINT(WMAX(J)+1D0)
        ELSE IF (FTYPE(J).EQ.'Y') THEN
*         Calendar axis.
          IF (GSTEP(J).LT.100D0) THEN
*           Daily increments.
            WMIN(J) = AINT(WMIN(J))
            WMAX(J) = AINT(WMAX(J)+1D0)
          ELSE
*           Start on Jan/01.
            CALL PGMJD (0, WMIN(J), IY, IM, ID)
            CALL PGMJD (1, WMIN(J), IY, 1, 1)

            CALL PGMJD (0, WMAX(J), IY, IM, ID)
            IF (GSTEP(J).LT.10000D0) THEN
*             Monthly increments.
              CALL PGMJD (1, WMAX(J), IY, 12, 1)
            ELSE
*             Yearly increments.
              CALL PGMJD (1, WMAX(J), IY+1, 1, 1)
            END IF
          END IF

        ELSE
          TMP = AINT(WMIN(J)/GSTEP(J))*GSTEP(J)
          IF (TMP.GE.WMIN(J)) TMP = TMP - GSTEP(J)
          WMIN(J) = TMP
          TMP = AINT(WMAX(J)/GSTEP(J))*GSTEP(J)
          IF (TMP.LE.WMAX(J)) TMP = TMP + GSTEP(J)
          WMAX(J) = TMP
        END IF

        DW(J) = WMAX(J) - WMIN(J)
        SW(J) = DW(J)/NSTEP(J)

*       Adjust NSTEP so that SW divides GSTEP.
        IF (SW(J).LT.GSTEP(J)) THEN
          SW(J) = GSTEP(J)/ANINT(GSTEP(J)/SW(J))
        ELSE
          SW(J) = GSTEP(J)
        END IF
        NSTEP(J) = ANINT(DW(J)/SW(J))
 70   CONTINUE


*  Draw the grid.
*     Get absolute scale for tick marks.
      CALL PGQVP (2, X1, X2, Y1, Y2)
      XSCL = (X2-X1)/(TRC(1)-BLC(1))
      YSCL = (Y2-Y1)/(TRC(2)-BLC(2))

*     Decode tick mark control.
      DO 80 J = 1, 2
        IF (GCODE(J).EQ.2) THEN
          TCODE(J,1) = -1
          TCODE(J,2) = -1
          TCODE(J,3) = -1
          TCODE(J,4) = -1
        ELSE IF (GCODE(J).EQ.1) THEN
          TCODE(J,1) = 1
          TCODE(J,2) = 1
          TCODE(J,3) = 1
          TCODE(J,4) = 1
        ELSE IF (GCODE(J).LT.0) THEN
          K = ABS(GCODE(J))
          TCODE(J,1) = MOD(K,10)
          TCODE(J,2) = MOD(K/10,10)
          TCODE(J,3) = MOD(K/100,10)
          TCODE(J,4) = MOD(K/1000,10)
        ELSE
          TCODE(J,1) = 0
          TCODE(J,2) = 0
          TCODE(J,3) = 0
          TCODE(J,4) = 0
        END IF
 80   CONTINUE

*     Draw each set of grid lines.
      OVERFL = .FALSE.
      DO 120 J = 1, 2
        IF (GCODE(J).EQ.0) GO TO 120

        IF (J.EQ.1) THEN
          CALL PGSCI (CJ(1))
          K = 2
        ELSE
          CALL PGSCI (CJ(2))
          K = 1
        END IF

        IF (NG(J).GT.0) THEN
          NWJ = NG(J)

        ELSE IF (FTYPE(J).EQ.'Y' .AND. GSTEP(J).GE.100D0) THEN
*         Calendar axis.
          CALL PGMJD (0, WMAX(J), IY, IM, ID)
          IF (GSTEP(J).LT.10000D0) THEN
            NWJ = 12*IY + IM
            CALL PGMJD (0, WMIN(J), IY, IM, ID)
            NWJ = (NWJ - (12*IY + IM))/INT(GSTEP(J)/100D0)
          ELSE
            NWJ = IY
            CALL PGMJD (0, WMIN(J), IY, IM, ID)
            NWJ = (NWJ - IY)/INT(GSTEP(J)/10000D0)
          END IF

        ELSE
          NWJ = NINT(DW(J)/GSTEP(J))
          IW0 = NINT(WMIN(J)/GSTEP(J))
        END IF

        DO 110 IWJ = 0, NWJ
          MAJOR = .FALSE.

*         Determine the world coordinate of the grid line.
          IF (NG(J).GT.0) THEN
*           User-specified.
            IF (IWJ.EQ.0) GO TO 110
            IF (J.EQ.1) THEN
              WORLD(1) = GRID1(IWJ)
            ELSE
              WORLD(2) = GRID2(IWJ)
            END IF

          ELSE
*           Internally computed.
            IF (FTYPE(J).EQ.'Y' .AND. GSTEP(J).GE.100D0) THEN
*             Calendar axis.
              CALL PGMJD (0, WMIN(J), IY, IM, ID)
              IF (GSTEP(J).LT.10000D0) THEN
                IM = IM + IWJ*INT(GSTEP(J)/100D0)
                CALL PGMJD (1, WORLD(J), IY, IM, ID)
              ELSE
                IY = IY + IWJ*INT(GSTEP(J)/10000D0)
                CALL PGMJD (1, WORLD(J), IY, IM, ID)
              END IF

            ELSE
              WORLD(J) = (IW0 + IWJ)*GSTEP(J)

*             Logarithmic?
              IF (FTYPE(J).EQ.'L') THEN
                TMP = MOD(WORLD(J),1D0)
                IF (TMP.LT.0D0) TMP = TMP + 1D0
                L = NINT(TMP*LDIV(J))

                IF (L.EQ.0) THEN
*                 Major tick mark.
                  MAJOR = .TRUE.
                ELSE
*                 Adjust logarithmic scales.
                  IF (LDIV(J).LE.6) THEN
                    L = LTABL(L,LDIV(J))
                  ELSE
                    L = L + 1
                  END IF

                  WORLD(J) = WORLD(J) - TMP + LOG10(DBLE(L))
                END IF
              END IF
            END IF
          END IF

          NP = 0
          GETEND = .TRUE.
          DO 100 IWK = 0, NSTEP(K)
            WORLD(K) = WMIN(K) + IWK*SW(K)

            IF (GETEND) THEN
*             Get end-point context.
              IF (DOLBOX) THEN
                XY(1) = (WORLD(1) - WXY(1))/XSPAN
                XY(2) = (WORLD(2) - WXY(3))/YSPAN
              ELSE
                CALL NLFUNC (1, NLC, NLI, NLD, NLCPRM, NLIPRM,
     :            NLDPRM, WORLD, XY, CONTRL, CONTXT, IERR)
              END IF

              IF (IERR.LE.0) THEN
                X1 = REAL(XY(1))
                Y1 = REAL(XY(2))
                INSIDE = X1.GT.BLC(1) .AND. X1.LT.TRC(1) .AND.
     :                   Y1.GT.BLC(2) .AND. Y1.LT.TRC(2)

                NP = 1
                XR(1) = X1
                YR(1) = Y1

                PREVIN = INSIDE
                GETEND = .FALSE.

                LABLOK = IERR.NE.-J .AND. IERR.NE.-3
              END IF
              GO TO 100
            END IF

            DO 90 ISTEP = 1, 1000
              IF (DOLBOX) THEN
                XY(1) = (WORLD(1) - WXY(1))/XSPAN
                XY(2) = (WORLD(2) - WXY(3))/YSPAN
              ELSE
                CALL NLFUNC (2, NLC, NLI, NLD, NLCPRM, NLIPRM,
     :            NLDPRM, WORLD, XY, CONTRL, CONTXT, IERR)
              END IF

              IF (IERR.GT.0) THEN
*               Flush buffer.
                IF (NP.GT.1) CALL PGLINE(NP, XR, YR)
                NP = 0
                GETEND = .TRUE.
                GO TO 100
              END IF

              IF (NP.EQ.BUFSIZ) THEN
*               Recycle buffer.
                CALL PGLINE(NP, XR, YR)
                XR(1) = XR(NP)
                YR(1) = YR(NP)
                NP = 1
              END IF

              X2 = REAL(XY(1))
              Y2 = REAL(XY(2))
              INSIDE = X2.GT.BLC(1) .AND. X2.LT.TRC(1) .AND.
     :                 Y2.GT.BLC(2) .AND. Y2.LT.TRC(2)

              IF (.NOT.INSIDE) THEN
*               For tick marks at the left or right edge.
                IF ((X2.EQ.BLC(1) .OR.  X2.EQ.TRC(1)) .AND.
     :               Y2.GT.BLC(2) .AND. Y2.LT.TRC(2)) THEN
                  INSIDE = X2.EQ.XR(NP)
                END IF
              END IF

              IF (.NOT.INSIDE) THEN
*               For tick marks at the bottom or top edge.
                IF ((Y2.EQ.BLC(2) .OR.  Y2.EQ.TRC(2)) .AND.
     :               X2.GT.BLC(1) .AND. X2.LT.TRC(1)) THEN
                  INSIDE = Y2.EQ.YR(NP)
                END IF
              END IF

              IF (NP.EQ.0) THEN
                NP = 1
                XR(1) = X2
                YR(1) = Y2
              ELSE
                IF (INSIDE) THEN
*                 This point is inside the frame...
                  IF (.NOT.PREVIN) THEN
*                   ...but the previous one was outside.
                    X1 = XR(NP)
                    Y1 = YR(NP)

                    FSEG = 0
                    XPOINT = 0.0
                    IF (ABS(X2-X1).GT.XTOL) THEN
                      S = (Y2-Y1)/(X2-X1)
                      IF (XR(NP).LE.BLC(1)) THEN
                        FSEG = 2
                        XR(NP) = BLC(1)
                        XPOINT = Y1 + (XR(NP) - X1)*S
                        YR(NP) = XPOINT
                      ELSE IF (XR(NP).GE.TRC(1)) THEN
                        FSEG = 4
                        XR(NP) = TRC(1)
                        XPOINT = Y1 + (XR(NP) - X1)*S
                        YR(NP) = XPOINT
                      END IF
                    END IF

                    IF (ABS(Y2-Y1).GT.YTOL) THEN
                      S = (X2-X1)/(Y2-Y1)
                      IF (YR(NP).LE.BLC(2)) THEN
                        FSEG = 1
                        YR(NP) = BLC(2)
                        XPOINT = X1 + (YR(NP) - Y1)*S
                        XR(NP) = XPOINT
                      ELSE IF (YR(NP).GE.TRC(2)) THEN
                        FSEG = 3
                        YR(NP) = TRC(2)
                        XPOINT = X1 + (YR(NP) - Y1)*S
                        XR(NP) = XPOINT
                      END IF
                    END IF


                    IF (FSEG.EQ.0) THEN
*                     The crossing is too oblique.
                      INSIDE = .FALSE.

                    ELSE
*                     Record this crossing point.
                      IF (TCODE(J,FSEG).NE.0 .AND. LABLOK) THEN
                        IF (IC.LT.NC-1) THEN
                          IC = IC + 1
                          CACHE(1,IC) = FSEG
                          CACHE(2,IC) = XPOINT
                          CACHE(3,IC) = J
                          CACHE(4,IC) = WORLD(J)
                        ELSE
*                         Cache overflow.
                          OVERFL = .TRUE.
                        END IF
                      END IF

                      IF (TCODE(J,FSEG).GT.0) THEN
*                       Just want tick marks.
                        S = (XSCL*(X2-X1))**2 + (YSCL*(Y2-Y1))**2
                        S = SQRT(S)/TCODE(J,FSEG)
                        IF (MAJOR) S = S/1.5
                        NP = NP + 1
                        XR(NP) = XR(NP-1) + (X2-X1)*TIKLEN/S
                        YR(NP) = YR(NP-1) + (Y2-Y1)*TIKLEN/S

                        CALL PGLINE(NP, XR, YR)
                        NP = 1
                      END IF
                    END IF
                  END IF

                  IF (INSIDE .AND. GCODE(J).EQ.2) THEN
*                   Full grid.
                    NP = NP + 1
                  END IF
                  XR(NP) = X2
                  YR(NP) = Y2
                ELSE
*                 This point is outside the frame...
                  IF (PREVIN) THEN
*                   ...but the previous one was inside.
                    X1 = XR(NP)
                    Y1 = YR(NP)

                    NP = NP + 1
                    XR(NP) = X2
                    YR(NP) = Y2

                    FSEG = 0
                    XPOINT = 0.0
                    IF (ABS(X2-X1).GT.XTOL) THEN
                      S = (Y2-Y1)/(X2-X1)
                      IF (XR(NP).LE.BLC(1)) THEN
                        FSEG = 2
                        XR(NP) = BLC(1)
                        XPOINT = Y1 + (XR(NP) - X1)*S
                        YR(NP) = XPOINT
                      ELSE IF (XR(NP).GE.TRC(1)) THEN
                        FSEG = 4
                        XR(NP) = TRC(1)
                        XPOINT = Y1 + (XR(NP) - X1)*S
                        YR(NP) = XPOINT
                      END IF
                    END IF

                    IF (ABS(Y2-Y1).GT.YTOL) THEN
                      S = (X2-X1)/(Y2-Y1)
                      IF (YR(NP).LE.BLC(2)) THEN
                        FSEG = 1
                        YR(NP) = BLC(2)
                        XPOINT = X1 + (YR(NP) - Y1)*S
                        XR(NP) = XPOINT
                      ELSE IF (YR(NP).GE.TRC(2)) THEN
                        FSEG = 3
                        YR(NP) = TRC(2)
                        XPOINT = X1 + (YR(NP) - Y1)*S
                        XR(NP) = XPOINT
                      END IF
                    END IF


                    IF (FSEG.EQ.0) THEN
*                     The crossing is too oblique.
                      INSIDE = .TRUE.

                      IF (GCODE(J).EQ.2) THEN
*                       Full grid.
                        NP = NP + 1
                      END IF
                      XR(NP) = X2
                      YR(NP) = Y2

                    ELSE
*                     Record this crossing point.
                      IF (TCODE(J,FSEG).NE.0 .AND. LABLOK) THEN
                        IF (IC.LT.NC-1) THEN
                          IC = IC + 1
                          CACHE(1,IC) = FSEG
                          CACHE(2,IC) = XPOINT
                          CACHE(3,IC) = J
                          CACHE(4,IC) = WORLD(J)
                        ELSE
*                         Cache overflow.
                          OVERFL = .TRUE.
                        END IF
                      END IF

                      IF (TCODE(J,FSEG).GT.0) THEN
*                       Just want tick marks.
                        X1 = XR(NP)
                        Y1 = YR(NP)
                        X2 = XR(NP-1)
                        Y2 = YR(NP-1)
                        S = (XSCL*(X2-X1))**2 + (YSCL*(Y2-Y1))**2
                        S = SQRT(S)/TCODE(J,FSEG)
                        IF (MAJOR) S = S/1.5
                        XR(NP-1) = X1 + (X2-X1)*TIKLEN/S
                        YR(NP-1) = Y1 + (Y2-Y1)*TIKLEN/S
                      END IF

*                     Flush buffer.
                      IF (TCODE(J,FSEG).NE.0) THEN
                        CALL PGLINE(NP, XR, YR)
                      END IF
                      NP = 0
                    END IF
                  ELSE
*                   The previous point was also outside.
                    XR(NP) = X2
                    YR(NP) = Y2
                  END IF
                END IF
              END IF

              PREVIN = INSIDE
              LABLOK = IERR.NE.-J .AND. IERR.NE.-3

              IF (CONTRL.EQ.0) THEN
                GO TO 100
              ELSE IF (CONTRL.EQ.1) THEN
*               Flush buffer.
                IF (NP.GT.1) CALL PGLINE(NP, XR, YR)
                NP = 0
              END IF
 90         CONTINUE
 100      CONTINUE

          IF (NP.GT.1) CALL PGLINE(NP, XR, YR)
 110    CONTINUE
 120  CONTINUE

      IERR = 0
      IF (OVERFL) IERR = 3


*  Produce axis labels.
 130  IF (LABCTL.NE.-1) CALL PGCRLB (BLC, TRC, IDENTS, FTYPE, LABCTL,
     :  CJ, NC, IC, CACHE)

*     Restore the original viewport, window and pen colour.
 999  CALL PGSVP (XVP1, XVP2, YVP1, YVP2)
      CALL PGSWIN (WXY(1), WXY(2), WXY(3), WXY(4))
      CALL PGSCI (CI0)

      RETURN
      END



*=======================================================================
*
* PGCRLB is a helper routine for PGSBOX, not meant to be called directly
* since it expects the viewport and window to be scaled to the full
* extent; it labels a curvilinear coordinate grid.
*
* Given:
*   BLC       R(2)      Cartesian coordinates of the bottom left-hand
*                       corner.
*   TRC       R(2)      Cartesian coordinates of the top right-hand
*                       corner.
*
*   IDENTS    C(3)*(*)  Identification strings (see PGSBOX).
*
*   FTYPE     C(2)*1    Axis types, used for axis labelling (see
*                       PGSBOX).
*
*   LABCTL    I         Decimal encoded grid labelling control (see
*                       PGSBOX).
*
*   CI        I(7)      Colour table (see PGSBOX).
*
* Given and/or returned:
*   NC        I         Upper array index for CACHE.
*
*   IC        I         Current number of entries in the CACHE table.
*
*   CACHE     D(4,0:NC) Table of points where the tick marks or grid
*                       lines cross the frame (see PGSBOX).
*
* Author: Mark Calabretta, Australia Telescope National Facility
*=======================================================================
      SUBROUTINE PGCRLB (BLC, TRC, IDENTS, FTYPE, LABCTL, CI, NC, IC,
     :                   CACHE)
*-----------------------------------------------------------------------
      LOGICAL   ANGLE, DODEG, DOMIN, DOYEAR, LFORCE, SEXA(2), TICKIT
      INTEGER   CI(7), DOLAB(4), EDGE, EDJE, IC, ID, ID2, IM, IM2,
     :          IMAG(2), ITER, IWRLD, IY, IY2, J, JC, K, K1, K2, KWRLD,
     :          L, LABCTL, LD, LM, LMAG(2), LS, LV, M, M1, M2, MM, NC,
     :          NCH, NI(2,0:4), NLAB, NLABS(2,-2:6), NSWAP, PP,
     :          PRVDEG(2), PRVMIN(2), PRVEDG, PRVYR(2), SEXSUP(2),
     :          SKOP(4)
      REAL      ANGL, BLC(2), FJUST, BNDRY(4), OMAG(2), SI(2), TRC(2),
     :          X, XBOX(4), XCH, XL, XW1, XW2, Y, YCH, YBOX(4), YL, YW1,
     :          YW2, Z
      DOUBLE PRECISION CACHE(4,0:NC), MJD1(2), MJD2(2), TMP, VS
      CHARACTER ESCAPE*1, EXPONT*20, FMT*8, IDENTS(3)*(*), TEXT*80,
     :          FTYPE(2)*1, TXT(2)*80

      DATA ESCAPE /'\\'/

*     These are to stop compiler messages about uninitialized variables.
      DATA DODEG, DOMIN, DOYEAR /3 * .FALSE./
      DATA XL, YL /2 * -999.0/
*-----------------------------------------------------------------------
*  Normalize angular table entries.
      IF (INDEX('ABDEGH',FTYPE(1)).NE.0 .OR.
     :    INDEX('ABDEGH',FTYPE(2)).NE.0) THEN
        DO 10 J = 1, IC
          IWRLD = NINT(CACHE(3,J))

          IF (INDEX('ADG', FTYPE(IWRLD)).NE.0) THEN
            CACHE(4,J) = MOD(CACHE(4,J), 360D0)
            IF (CACHE(4,J).LT.0D0) CACHE(4,J) = CACHE(4,J) + 360D0
          ELSE IF (INDEX('BEH', FTYPE(IWRLD)).NE.0) THEN
            CACHE(4,J) = MOD(CACHE(4,J), 360D0)
            IF (CACHE(4,J).LE.-180D0) THEN
              CACHE(4,J) = CACHE(4,J) + 360D0
            ELSE IF (CACHE(4,J).GT.180D0) THEN
              CACHE(4,J) = CACHE(4,J) - 360D0
            END IF
          END IF

          IF (INDEX('GHI', FTYPE(IWRLD)).NE.0) THEN
*           Angle expressed as time.
            CACHE(4,J) = CACHE(4,J)/15D0
          END IF
 10     CONTINUE
      END IF


*  Reorganize the table entries.
*     Sort crossings for each of the four frame segments.
      DO 40 ITER = 1, IC
        NSWAP = 0
        DO 30 J = 1, IC-1
          IF (CACHE(1,J).LT.CACHE(1,J+1)) GO TO 30
          IF (CACHE(1,J).EQ.CACHE(1,J+1) .AND.
     :        CACHE(2,J).LE.CACHE(2,J+1)) GO TO 30

          NSWAP = NSWAP + 1
          DO 20 M = 1, 4
            TMP = CACHE(M,J)
            CACHE(M,J) = CACHE(M,J+1)
            CACHE(M,J+1) = TMP
 20       CONTINUE
 30     CONTINUE
        IF (NSWAP.EQ.0) GO TO 50
 40   CONTINUE

*     Squeeze out duplicates.
 50   JC = IC
      DO 90 J = 2, IC
        IF (J.GT.JC) GO TO 100

        DO 60 M = 1, 4
          IF (CACHE(M,J).NE.CACHE(M,J-1)) GO TO 90
 60     CONTINUE

*       This entry is the same as the previous one.
        JC = JC - 1
        DO 80 K = J, JC
          DO 70 M = 1, 4
            CACHE(M,K) = CACHE(M,K+1)
 70       CONTINUE
 80     CONTINUE
 90   CONTINUE

 100  IC = JC


* How do we label the edges of the frame?
*     Determine separability indices.
      NI(1,0) = 0
      NI(1,1) = 0
      NI(1,2) = 0
      NI(1,3) = 0
      NI(1,4) = 0
      NI(2,0) = 0
      NI(2,1) = 0
      NI(2,2) = 0
      NI(2,3) = 0
      NI(2,4) = 0
      DO 110 J = 1, IC
        IWRLD = NINT(CACHE(3,J))
        EDGE  = NINT(CACHE(1,J))

        NI(IWRLD,0) = NI(IWRLD,0) + 1
        NI(IWRLD,EDGE) = NI(IWRLD,EDGE) + 1
 110  CONTINUE

      SI(1) = 0.0
      SI(2) = 0.0
      IF (NI(1,0).GT.0) SI(1) = 2.0*REAL(NI(1,1)+NI(1,3))/NI(1,0) - 1.0
      IF (NI(2,0).GT.0) SI(2) = 2.0*REAL(NI(2,1)+NI(2,3))/NI(2,0) - 1.0

*     Which coordinates go on which edges?
      IF (LABCTL.GT.0) THEN
*       User-defined.
        L = LABCTL
      ELSE
*       Work it out ourselves.
        L = 0
        IF (ABS(SI(1)-SI(2)).GT.1.0) THEN
*         Approximately horizontal/vertical grid lines.
          IF (SI(1).GT.SI(2)) THEN
*           First world coordinate with vertical grid lines.
            IF (NI(1,1).GT.3 .OR. NI(1,1).GE.NI(1,3)) THEN
*             Label bottom of frame.
              L = L + 1
            ELSE
*             Label top of frame.
              L = L + 100
            END IF

*           Second world coordinate with horizontal grid lines.
            IF (NI(2,2).GT.3 .OR. NI(2,2).GE.NI(2,4)) THEN
*             Label left side of frame.
              L = L + 20
            ELSE
*             Label right side of frame.
              L = L + 2000
            END IF
          ELSE
*           First world coordinate with horizontal grid lines.
            IF (NI(1,2).GT.3 .OR. NI(1,2).GE.NI(1,4)) THEN
*             Label left side of frame.
              L = L + 10
            ELSE
*             Label right side of frame.
              L = L + 1000
            END IF

*           Second world coordinate with vertical grid lines.
            IF (NI(2,1).GT.3 .OR. NI(2,1).GE.NI(2,3)) THEN
*             Label bottom of frame.
              L = L + 2
            ELSE
*             Label top of frame.
              L = L + 200
            END IF
          END IF
        ELSE
*         Skew grid lines or worse.
          IF (SI(1).GT.0.5D0) THEN
*           First world coordinate with vertical grid lines.
            IF (NI(1,1).GT.3 .OR. NI(1,1).GT.NI(1,3)) THEN
              L = 1
            ELSE
              L = 100
            END IF

            IF (SI(2).GT.0.5D0) THEN
*             Second world coordinate also with vertical grid lines.
              IF (L.EQ.1) THEN
                L = 201
              ELSE
                IF (NI(2,1).GT.1 .OR. NI(2,1).GT.NI(2,3)) THEN
                  L = 102
                ELSE
                  L = 300
                END IF
              END IF
            ELSE
*             Second world coordinate with diagonal grid lines.
              L = L + 2020
            END IF
          ELSE IF (SI(1).LT.-0.5D0) THEN
*           First world coordinate with horizontal grid lines.
            IF (NI(1,2).GT.3 .OR. NI(1,2).GT.NI(1,4)) THEN
              L = 10
            ELSE
              L = 1000
            END IF

            IF (SI(2).LT.-0.5D0) THEN
*             Second world coordinate also with horizontal grid.
              IF (L.EQ.10) THEN
                L = 2010
              ELSE
                IF (NI(2,2).GT.1 .OR. NI(2,2).GT.NI(2,4)) THEN
                  L = 1020
                ELSE
                  L = 3000
                END IF
              END IF
            ELSE
*             Second world coordinate with diagonal grid lines.
              L = L + 202
            END IF
          ELSE IF (SI(2).GT.0.5D0) THEN
*           Second world coordinate with vertical grid lines.
            IF (NI(2,1).GT.3 .OR. NI(2,1).GT.NI(2,3)) THEN
              L = 1012
            ELSE
              L = 1210
            END IF
          ELSE IF (SI(2).LT.-0.5D0) THEN
*           Second world coordinate with horizontal grid lines.
            IF (NI(2,2).GT.3 .OR. NI(2,2).GT.NI(2,4)) THEN
              L = 121
            ELSE
              L = 2101
            END IF
          ELSE
*           Desperation stakes!  Label all four axes.
            K1 = NI(1,3) + NI(1,1)
            K2 = NI(2,4) + NI(2,2)
            L = 2121

            M1 = NI(1,4) + NI(1,1)
            M2 = NI(2,3) + NI(2,2)
            IF (M1.GE.K1 .AND. M2.GE.K2) THEN
              L  = 1221
              K1 = M1
              K2 = M2
            END IF

            M1 = NI(1,2) + NI(1,1)
            M2 = NI(2,4) + NI(2,3)
            IF (M1.GE.K1 .AND. M2.GE.K2) THEN
              L  = 2211
              K1 = M1
              K2 = M2
            END IF

            M1 = NI(1,4) + NI(1,2)
            M2 = NI(2,3) + NI(2,1)
            IF (M1.GE.K1 .AND. M2.GE.K2) THEN
              L  = 1212
              K1 = M1
              K2 = M2
            END IF

            M1 = NI(1,3) + NI(1,2)
            M2 = NI(2,4) + NI(2,1)
            IF (M1.GE.K1 .AND. M2.GE.K2) THEN
              L  = 2112
              K1 = M1
              K2 = M2
            END IF

            M1 = NI(1,4) + NI(1,3)
            M2 = NI(2,2) + NI(2,1)
            IF (M1.GE.K1 .AND. M2.GE.K2) THEN
              L  = 1122
              K1 = M1
              K2 = M2
            END IF
          END IF
        END IF
      END IF

*     Mark labels as unwanted by setting their edges negative.
      DOLAB(1) = MOD(L,10)
      DOLAB(2) = MOD(L/10,10)
      DOLAB(3) = MOD(L/100,10)
      DOLAB(4) = MOD(L/1000,10)
      DO 120 J = 1, IC
        EDGE  = NINT(CACHE(1,J))
        IWRLD = NINT(CACHE(3,J))
        IF (IWRLD.EQ.1) THEN
          IF (MOD(DOLAB(EDGE),2).NE.1) CACHE(1,J) = -EDGE
        ELSE
          IF (MOD(DOLAB(EDGE)/2,2).NE.1) CACHE(1,J) = -EDGE
        END IF
 120  CONTINUE


*  Determine labelling precision.
*     Order of magnitude for plain numeric world coordinates.
      IMAG(1) = 0
      IMAG(2) = 0

      IF (FTYPE(1).EQ.' ' .OR. FTYPE(2).EQ.' ') THEN
        OMAG(1) = 0.0
        OMAG(2) = 0.0
        DO 130 J = 1, IC
          IWRLD = NINT(CACHE(3,J))
          IF (FTYPE(IWRLD).EQ.' ') THEN
*           Plain numeric.
            IF (CACHE(4,J).EQ.0D0) GO TO 130
            OMAG(IWRLD) = OMAG(IWRLD) + LOG10(ABS(CACHE(4,J)))
            IMAG(IWRLD) = IMAG(IWRLD) + 1
          END IF
 130    CONTINUE

        IF (IMAG(1).GT.0) IMAG(1) = INT(OMAG(1)/IMAG(1))
        IF (IMAG(1).GE.-2 .AND. IMAG(1).LE.4) IMAG(1) = 0

        IF (IMAG(2).GT.0) IMAG(2) = INT(OMAG(2)/IMAG(2))
        IF (IMAG(2).GE.-2 .AND. IMAG(2).LE.4) IMAG(2) = 0

*       Renormalize grid values.
        IF (IMAG(1).NE.0 .OR. IMAG(2).NE.0) THEN
          DO 140 J = 1, IC
            IWRLD = NINT(CACHE(3,J))
            CACHE(4,J) = CACHE(4,J)/10D0**IMAG(IWRLD)
 140      CONTINUE
        END IF
      END IF

*     Sexagesimal labelling.
      SEXA(1) = INDEX('DEFGHITy', FTYPE(1)).NE.0
      SEXA(2) = INDEX('DEFGHITy', FTYPE(2)).NE.0

      IF (SEXA(1) .OR. SEXA(2)) THEN
        DO 150 K = -2, 6
          NLABS(1,K) = 0
          NLABS(2,K) = 0
 150    CONTINUE

        DO 180 J = 1, IC
*         Use this label?
          EDGE = NINT(CACHE(1,J))
          IF (EDGE.LT.0) GO TO 180

*         Skip non-sexagesimal coordinates.
          IWRLD = NINT(CACHE(3,J))
          IF (.NOT.SEXA(IWRLD)) GO TO 180

*         Defeat rounding errors.
          IF (FTYPE(IWRLD).EQ.'y') THEN
            TMP = ABS(MOD(CACHE(4,J),1D0)*86400D0) + 5D-7
          ELSE
            TMP = ABS(CACHE(4,J)*3600D0) + 5D-7
          END IF

          LV = INT(MOD(TMP, 1D0)*1D6)
          IF (LV.NE.0) THEN
*           Sub-arcsec/second resolution.
            M = 1
            DO 160 K = 1, 6
              M = M*10
              IF (MOD(LV,M).NE.0) THEN
                L = 7 - K
                GO TO 170
              END IF
 160        CONTINUE
          ELSE
            LS = INT(TMP)
            IF (MOD(LS,60).NE.0) THEN
              L =  0
            ELSE IF (MOD(LS,3600).NE.0) THEN
              L = -1
            ELSE
              L = -2
            END IF
          END IF

*         Use CACHE(1,J) temporarily to hold L.
 170      NLABS(IWRLD,L) = NLABS(IWRLD,L) + 1
          CACHE(1,J) = CACHE(1,J) + 0.1 + L/100.0
 180    CONTINUE

*       How many fields are needed to get enough labels?
        DO 200 IWRLD = 1, 2
          NLAB = 0
          DO 190 K = -2, 6
            NLAB = NLAB + NLABS(IWRLD,K)
            IF (NLAB.GT.3 .OR. K.EQ.6) THEN
              LMAG(IWRLD) = K
              GO TO 200
            END IF
            NLABS(IWRLD,K+1) = NLABS(IWRLD,K) + NLABS(IWRLD,K+1)
 190      CONTINUE
 200    CONTINUE

*       Disable unwanted labels.
        DO 210 J = 1, IC
*         Use this label?
          EDGE = NINT(CACHE(1,J))
          IF (EDGE.LT.0) GO TO 210

*         Skip non-sexagesimal coordinates.
          IWRLD = NINT(CACHE(3,J))
          IF (.NOT.SEXA(IWRLD)) GO TO 210

          L = NINT(100*MOD(CACHE(1,J),1.0)) - 10
          IF (L.LE.LMAG(IWRLD)) THEN
*           Wanted.
            CACHE(1,J) = EDGE
          ELSE
*           Not wanted.
            CACHE(1,J) = -EDGE
          END IF
 210    CONTINUE
      END IF


*  Produce labels.
      XW1 = BLC(1)
      XW2 = TRC(1)
      YW1 = BLC(2)
      YW2 = TRC(2)

*     These will define a box that just contains the labels.
      BNDRY(1) = YW1
      BNDRY(2) = XW1
      BNDRY(3) = YW2
      BNDRY(4) = XW2

*     These will record the edges on which labels were actually written.
      SKOP(1) = 0
      SKOP(2) = 0
      SKOP(3) = 0
      SKOP(4) = 0

*     Calendar date range.
      MJD1(1) =  1D99
      MJD1(2) =  1D99
      MJD2(1) = -1D99
      MJD2(2) = -1D99

*     Character height.
      CALL PGQCS (4, XCH, YCH)
      PRVEDG = 0

*     Loop through the axis crossing table.
      DO 340 J = 1, IC
*       Determine the position.
        IWRLD = NINT(CACHE(3,J))
        CALL PGSCI (CI(IWRLD+2))

        EDGE = ABS(NINT(CACHE(1,J)))
        IF (EDGE.NE.PRVEDG) THEN
*         Start of new edge.

          IF (SEXA(1) .OR. SEXA(2)) THEN
*           Sexagesimal field suppression policy.
            PRVDEG(1) = -1
            PRVMIN(1) = -1
            PRVDEG(2) = -1
            PRVMIN(2) = -1

            SEXSUP(1) = 0
            SEXSUP(2) = 0

*           Vertical sides.
            IF (MOD(EDGE,2).EQ.0) THEN
              DO 220 K = J, IC
                EDJE = ABS(NINT(CACHE(1,K)))
                IF (EDJE.NE.EDGE) GO TO 230

                KWRLD = NINT(CACHE(3,K))
                IF (SEXSUP(KWRLD).EQ.1) GO TO 220

                IF (FTYPE(KWRLD).EQ.'y') THEN
                  TMP = ABS(MOD(CACHE(4,K),1D0)*24D0)
                ELSE
                  TMP = ABS(CACHE(4,K))
                END IF

                LV = INT(TMP*3600D0 + 5D-7)
                LD =  LV/3600
                LM = (LV - LD*3600)/60
                LS =  LV - LD*3600 - LM*60
                TMP = TMP - LD

                IF (TMP.LT.5D-7) THEN
*                 Write deg/hr field only when min/sec are zero.
                  SEXSUP(KWRLD) = 1
                ELSE IF (TMP*60-LM.LT.3D-5) THEN
*                 Write min field only when sec is zero; only
*                 write the deg/hr field when min is written.
                  SEXSUP(KWRLD) = 2
                END IF
 220          CONTINUE
            END IF
          END IF
 230      CONTINUE

          PRVYR(1) = -1
          PRVYR(2) = -1

          XL = -999.0
          YL = -999.0
        END IF
        PRVEDG = EDGE
        LFORCE = .FALSE.

*       Fix up edge encoding.
        EDGE = NINT(CACHE(1,J))
        CACHE(1,J) = ABS(EDGE)

*       Use this label?
        IF (EDGE.LT.0) GO TO 340

        IF (EDGE.EQ.1) THEN
*         Bottom.
          FJUST = 0.5
          X = CACHE(2,J)
          Y = YW1 - 1.5*YCH
        ELSE IF (EDGE.EQ.2) THEN
*         Left.
          FJUST = 1.0
          X = XW1 - 0.5*XCH
          Y = CACHE(2,J) - YCH/2.0
        ELSE IF (EDGE.EQ.3) THEN
*         Top.
          FJUST = 0.5
          X = CACHE(2,J)
          Y = YW2 + 0.5*YCH
        ELSE IF (EDGE.EQ.4) THEN
*         Right.
          FJUST = 0.0
          X = XW2 + 0.5*XCH
          Y = CACHE(2,J) - YCH/2.0
        END IF

*       Format the numeric label.
        IF (INDEX('ABC', FTYPE(IWRLD)).NE.0) THEN
*         Decimal angle; allow up to 6 decimal digits.
          TMP = ABS(CACHE(4,J)) + 5D-7
          LD  = INT(TMP)

          K = 1
          IF (CACHE(4,J).LT.0D0) THEN
*           Insert a minus sign.
            TEXT(1:1) = '-'
            K = 2
          END IF

          CALL PGNUMB (LD, 0, 1, TEXT(K:), NCH)
          K = K + NCH
          TEXT(K:) = ESCAPE // 'uo' // ESCAPE // 'd'
          K = K + 4

          LV = INT(MOD(TMP,1D0)*1D6)
          IF (LV.NE.0) CALL PGNUMB (LV, -6, 1, TEXT(K:), NCH)
          TEXT(K:K) = 'd'

        ELSE IF (SEXA(IWRLD)) THEN
*         Sexagesimal format; angle or time?
          ANGLE = INDEX('DEF', FTYPE(IWRLD)).NE.0
          L = LMAG(IWRLD)

*         Use integer arithmetic to avoid rounding problems.
          IF (FTYPE(IWRLD).EQ.'y') THEN
            TMP = ABS(MOD(CACHE(4,J),1D0)*24D0)

*           Determine date range.
            IF (CACHE(4,J).LT.MJD1(IWRLD)) MJD1(IWRLD) = CACHE(4,J)
            IF (CACHE(4,J).GT.MJD2(IWRLD)) MJD2(IWRLD) = CACHE(4,J)
          ELSE
            TMP = ABS(CACHE(4,J))
          END IF
          VS = TMP*3600D0 + 5D-7
          LV = INT(VS)

*         Sexagesimal fields.
          LD =  LV/3600
          LM = (LV - LD*3600)/60
          LS =  LV - LD*3600 - LM*60

*         Field suppression policy.
          IF (SEXSUP(IWRLD).GT.0) THEN
            TMP = TMP - LD

            IF (TMP.LT.5D-7) THEN
              DODEG = .TRUE.
              DOMIN = .TRUE.
            ELSE IF (SEXSUP(IWRLD).EQ.2) THEN
              DOMIN = TMP*60-LM.LT.3D-5
              DODEG = DOMIN .AND. LD.NE.PRVDEG(IWRLD)
            ELSE
              DODEG = .FALSE.
              DOMIN = .FALSE.
            END IF
          ELSE
            DODEG = LD.NE.PRVDEG(IWRLD)
            DOMIN = LM.NE.PRVMIN(IWRLD)
          END IF

          K = 1
          IF (L.EQ.-2 .OR. DODEG) THEN
*           Write the degree/hour field.
            DODEG = .TRUE.

            IF (CACHE(4,J).LT.0D0) THEN
*             Insert a minus sign.
              TEXT(1:1) = '-'
              K = 2
            END IF

            IF (ANGLE) THEN
*             Angle.
              CALL PGNUMB (LD, 0, 1, TEXT(K:), NCH)
              K = K + NCH
              TEXT(K:) = ESCAPE // 'uo' // ESCAPE // 'd'
              K = K + 5
            ELSE
*             Time.
              IF (LD.LE.9 .AND. INDEX('Ty',FTYPE(IWRLD)).EQ.0) THEN
*               Write leading zeroes in the hour field.
                WRITE (TEXT(K:), '(I2.2)') LD
                K = K + 2
              ELSE
                CALL PGNUMB (LD, 0, 1, TEXT(K:), NCH)
                K = K + NCH
              END IF
              TEXT(K:) = ESCAPE // 'uh' // ESCAPE // 'd'
              K = K + 5
            END IF
          END IF

          IF (L.GE.-1 .AND.
     :      .NOT.(L.LE.0 .AND. K.GT.1 .AND. LM.EQ.0 .AND. LS.EQ.0)) THEN
*           Write arcminute/minute field.
            IF (L.EQ.-1 .OR. K.GT.1 .OR. DOMIN) THEN
              DOMIN = .TRUE.
              IF (ANGLE) THEN
                WRITE (TEXT(K:), '(I2.2,A)') LM, ''''
                K = K + 3
              ELSE
                WRITE (TEXT(K:), '(I2.2,A)') LM,
     :            ESCAPE // 'um' // ESCAPE // 'd'
                K = K + 7
              END IF
            END IF

            IF (L.GE.0) THEN
*             Arcsec/second field.
              IF (ANGLE) THEN
                WRITE (TEXT(K:), '(I2.2,A)') LS, '"'
                K = K + 3
              ELSE
                WRITE (TEXT(K:), '(I2.2,A)') LS,
     :            ESCAPE // 'us' // ESCAPE // 'd'
                K = K + 7
              END IF

              IF (L.GT.0) THEN
*               Sub-arcsec/second field.
                WRITE (FMT, '(A,I1,A,I1,A)') '(A,I', L, '.', L, ')'
                LV = INT(MOD(VS,1D0)*10**L)
                WRITE (TEXT(K:), FMT) '.', LV
              END IF
            END IF
          END IF

        ELSE IF (FTYPE(IWRLD).EQ.'L') THEN
*         Logarithmic.
          TMP = MOD(CACHE(4,J),1D0)
          IF (TMP.LT.0D0) TMP = TMP + 1D0
          MM = NINT(10D0**TMP)
          IF (MM.EQ.10) THEN
            TMP = TMP - 1D0
            MM = 1
          END IF

          PP = NINT(CACHE(4,J) - TMP)

          IF (MM.NE.1) THEN
            WRITE (TEXT, '(I1)') MM
          ELSE
*           FORTRAN is really abysmal sometimes.
            WRITE (TEXT, '(I8)') PP
            DO 240 K = 1, 8
              IF (TEXT(K:K).NE.' ') GO TO 250
 240        CONTINUE
 250        TEXT = '10' // ESCAPE // 'u' // TEXT(K:8)

            LFORCE = .TRUE.
          END IF

        ELSE IF (FTYPE(IWRLD).EQ.'Y') THEN
*         Convert MJD to Gregorian calendar date, YYYY/MM/DD.
          CALL PGMJD(0, CACHE(4,J), IY, IM, ID)
          DOYEAR = IY.NE.PRVYR(IWRLD)

          IF (DOYEAR) THEN
            WRITE (TEXT, 260) IY, IM, ID
 260        FORMAT (I12,'/',I2.2,'/',I2.2)

            DO 270 K = 1, 12
              IF (TEXT(K:K).NE.' ') GO TO 280
 270        CONTINUE
 280        TEXT = TEXT(K:18)

          ELSE
            WRITE (TEXT, 290) IM, ID
 290        FORMAT (I2.2,'/',I2.2)
          END IF

        ELSE
*         Plain number; allow up to six significant digits.
          IF (CACHE(4,J).NE.0D0) THEN
            PP = INT(LOG10(ABS(CACHE(4,J)))) - 6
            MM = NINT(CACHE(4,J)/10D0**PP)
          ELSE
            PP = 0
            MM = 0
          END IF

          CALL PGNUMB (MM, PP, 0, TEXT, NCH)
        END IF

*       Write the label if it doesn't overlap the previous one.
        CALL PGQTXT (X, Y, 0.0, FJUST, TEXT, XBOX, YBOX)
        IF (LFORCE .OR. XBOX(1).GT.XL .OR. YBOX(1).GT.YL) THEN
          IF (IWRLD.EQ.1) THEN
            CALL PGSCI (CI(3))
          ELSE
            CALL PGSCI (CI(4))
          END IF

          CALL PGPTXT (X, Y, 0.0, FJUST, TEXT)
          XL = XBOX(4) + 0.5*XCH
          YL = YBOX(2) + 0.5*YCH

*         Sexagesimal formatting.
          IF (SEXA(IWRLD)) THEN
            IF (DODEG) PRVDEG(IWRLD) = LD
            IF (DOMIN) PRVMIN(IWRLD) = LM
          END IF

*         Calendar formatting.
          IF (FTYPE(IWRLD).EQ.'Y') THEN
            IF (DOYEAR) PRVYR(IWRLD) = IY
          END IF

*         Record the fact.
          IF (IWRLD.EQ.1) THEN
            SKOP(EDGE) = 2*(SKOP(EDGE)/2) + 1
          ELSE
            SKOP(EDGE) = 2 + MOD(SKOP(EDGE),2)
          END IF

*         Boundary within which the numeric labels lie.
          IF (YBOX(1).LT.BNDRY(1)) BNDRY(1) = YBOX(1)
          IF (XBOX(1).LT.BNDRY(2)) BNDRY(2) = XBOX(1)
          IF (YBOX(3).GT.BNDRY(3)) BNDRY(3) = YBOX(3)
          IF (XBOX(3).GT.BNDRY(4)) BNDRY(4) = XBOX(3)

*         Check the distance to the previous grid line.
          TICKIT = .FALSE.
          K = J - 1
          IF (DOLAB(EDGE).NE.3) THEN
*           Only one coordinate element labelled on this edge, no need
*           to worry about grid lines belonging to the other.
            DO 300 K = J-1, 1, -1
              IF (IWRLD.EQ.NINT(CACHE(3,K))) GO TO 310
 300        CONTINUE
          END IF

 310      IF (K.GE.1) THEN
            EDJE = ABS(NINT(CACHE(1,K)))
            IF (EDJE.EQ.EDGE) THEN
              IF (EDGE.EQ.1 .OR. EDGE.EQ.3) THEN
                IF (XBOX(1).LT.CACHE(2,K)) TICKIT = .TRUE.
              ELSE
                IF (YBOX(1).LT.CACHE(2,K)) TICKIT = .TRUE.
              END IF
            END IF
          END IF

*         Check the distance to the next grid line.
          K = J + 1
          IF (DOLAB(EDGE).NE.3) THEN
*           Only one coordinate element labelled on this edge, no need
*           to worry about grid lines belonging to the other.
            DO 320 K = J+1, IC
              IF (IWRLD.EQ.NINT(CACHE(3,K))) GO TO 330
 320        CONTINUE
          END IF

 330      IF (K.LE.IC) THEN
            EDJE = ABS(NINT(CACHE(1,K)))
            IF (EDJE.EQ.EDGE) THEN
              IF (EDGE.EQ.1 .OR. EDGE.EQ.3) THEN
                IF (XBOX(3).GT.CACHE(2,K)) TICKIT = .TRUE.
              ELSE
                IF (YBOX(3).GT.CACHE(2,K)) TICKIT = .TRUE.
              END IF
            END IF
          END IF

*         Density of grid lines is high.
          IF (TICKIT) THEN
*           Draw outside pip mark to disambiguate grid lines.
            Z = REAL(CACHE(2,J))
            IF (EDGE.EQ.1) THEN
              CALL PGMOVE (Z, YW1-0.3*YCH)
              CALL PGDRAW (Z, YW1)
            ELSE IF (EDGE.EQ.2) THEN
              CALL PGMOVE (XW1-0.3*XCH, Z)
              CALL PGDRAW (XW1, Z)
            ELSE IF (EDGE.EQ.3) THEN
              CALL PGMOVE (Z, YW2+0.3*YCH)
              CALL PGDRAW (Z, YW2)
            ELSE IF (EDGE.EQ.4) THEN
              CALL PGMOVE (XW2+0.3*XCH, Z)
              CALL PGDRAW (XW2, Z)
            END IF
          END IF
        END IF
 340  CONTINUE


*  Write the identification strings.
*     World coordinates.
      DO 470 EDGE = 1, 4
        TEXT = ' '

        DO 440 IWRLD = 1, 2
          TXT(IWRLD) = ' '

          IF (MOD(SKOP(EDGE)/IWRLD,2).EQ.0) GO TO 440

*         Strip off leading blanks.
          L = LEN(IDENTS(IWRLD))
          DO 350 K = 1, L
            IF (IDENTS(IWRLD)(K:K).NE.' ') THEN
              TXT(IWRLD) = IDENTS(IWRLD)(K:L)
              GO TO 360
            END IF
 350      CONTINUE

 360      IF (IMAG(IWRLD).NE.0 .OR. FTYPE(IWRLD).EQ.'y') THEN
*           Find the last non-blank.
            DO 370 K = 40, 1, -1
              IF (TXT(IWRLD)(K:K).NE.' ') GO TO 380
 370        CONTINUE
 380        K = K + 1

            IF (IMAG(IWRLD).NE.0) THEN
*             Add scaling information.
              CALL PGNUMB (IMAG(IWRLD), 0, 1, EXPONT, NCH)
              TXT(IWRLD)(K:) = '  x10' // ESCAPE // 'u' // EXPONT
            ELSE
*             Add calendar date range.
              CALL PGMJD(0, MJD1(IWRLD), IY, IM, ID)

              WRITE (TEXT, 390) IY, IM, ID
 390          FORMAT (I12,'/',I2.2,'/',I2.2)
              DO 400 L = 1, 12
                IF (TEXT(L:L).NE.' ') GO TO 410
 400          CONTINUE
 410          TXT(IWRLD)(K:) = ' (' // TEXT(L:18)
              K = K + 21 - L

              CALL PGMJD(0, MJD2(IWRLD), IY2, IM2, ID2)
              WRITE (TEXT, 390) IY2, IM2, ID2
              IF (IY2.EQ.IY) THEN
                IF (IM2.EQ.IM) THEN
                  IF (ID2.EQ.ID) THEN
                    TXT(IWRLD)(K:) = ')'
                  ELSE
                    TXT(IWRLD)(K:) = ' - ' // TEXT(17:18) // ')'
                  END IF
                ELSE
                  TXT(IWRLD)(K:) = ' - ' // TEXT(14:18) // ')'
                END IF
              ELSE
                DO 420 L = 1, 12
                  IF (TEXT(L:L).NE.' ') GO TO 430
 420            CONTINUE
 430            TXT(IWRLD)(K:) = ' - ' // TEXT(L:18) // ')'
              END IF
            END IF
          END IF

 440    CONTINUE

        K = 0
        IF (TXT(1).NE.' ') THEN
*         Identify first world coordinate...
          TEXT = TXT(1)
          CALL PGSCI (CI(5))

          IF (TXT(2).NE.' ') THEN
*           ...and also second world coordinate.
            DO 450 K = 40, 1, -1
              IF (TEXT(K:K).NE.' ') GO TO 460
 450        CONTINUE
 460        K = K + 1

            TEXT(K:) = ',  ' // TXT(2)
          END IF
        ELSE IF (TXT(2).NE.' ') THEN
*         Identify second world coordinate only.
          TEXT = TXT(2)
          CALL PGSCI (CI(6))
        ELSE
*         No text to write.
          GO TO 470
        END IF

        IF (EDGE.EQ.1) THEN
          X = (XW1 + XW2)/2.0
          Y = BNDRY(1) - 1.5*YCH
          ANGL = 0.0
        ELSE IF (EDGE.EQ.2) THEN
          IF (TEXT(2:).EQ.' ') THEN
*           One character, write upright.
            X = BNDRY(2) - 1.0*XCH
            Y = (YW1 + YW2)/2.0
            ANGL = 0.0
          ELSE
            X = BNDRY(2) - 0.5*XCH
            Y = (YW1 + YW2)/2.0
            ANGL = 90.0
          END IF
        ELSE IF (EDGE.EQ.3) THEN
          X = (XW1 + XW2)/2.0
          Y = BNDRY(3) + 0.5*YCH
          ANGL = 0.0
        ELSE IF (EDGE.EQ.4) THEN
          IF (TEXT(2:).EQ.' ') THEN
*           One character, write upright.
            X = BNDRY(4) + 1.0*XCH
            Y = (YW1 + YW2)/2.0
            ANGL = 0.0
          ELSE
            X = BNDRY(4) + 0.5*XCH
            Y = (YW1 + YW2)/2.0
            ANGL = -90.0
          END IF
        END IF

        CALL PGQTXT (X, Y, ANGL, 0.5, TEXT, XBOX, YBOX)
        IF (K.EQ.0) THEN
          CALL PGPTXT (X, Y, ANGL, 0.5, TEXT)
        ELSE
*         Two-colour annotation.
          CALL PGPTXT (XBOX(1), YBOX(1), ANGL, 0.0, TEXT(:K))
          CALL PGSCI (CI(6))
          CALL PGPTXT (XBOX(4), YBOX(4), ANGL, 1.0, TEXT(K+1:))
        END IF

*       Update the boundary.
        IF (YBOX(1).LT.BNDRY(1)) BNDRY(1) = YBOX(1)
        IF (YBOX(3).LT.BNDRY(1)) BNDRY(1) = YBOX(3)
        IF (XBOX(1).LT.BNDRY(2)) BNDRY(2) = XBOX(1)
        IF (XBOX(3).LT.BNDRY(2)) BNDRY(2) = XBOX(3)
        IF (YBOX(1).GT.BNDRY(3)) BNDRY(3) = YBOX(1)
        IF (YBOX(3).GT.BNDRY(3)) BNDRY(3) = YBOX(3)
        IF (XBOX(1).GT.BNDRY(4)) BNDRY(4) = XBOX(1)
        IF (XBOX(3).GT.BNDRY(4)) BNDRY(4) = XBOX(3)
 470  CONTINUE

*     Title.
      IF (IDENTS(3).NE.' ') THEN
        CALL PGSCI (CI(7))
        X = (XW1 + XW2)/2.0
        Y = BNDRY(3) + 0.5*YCH
        CALL PGQTXT (X, Y, 0.0, 0.5, IDENTS(3), XBOX, YBOX)
        CALL PGPTXT (X, Y, 0.0, 0.5, IDENTS(3))

*       Update the boundary.
        IF (XBOX(1).LT.BNDRY(2)) BNDRY(2) = XBOX(1)
        IF (YBOX(3).GT.BNDRY(3)) BNDRY(3) = YBOX(3)
        IF (XBOX(3).GT.BNDRY(4)) BNDRY(4) = XBOX(3)
      END IF


*     Return margin widths in Cartesian coordinates.
      CACHE(1,NC) = YW1 - BNDRY(1)
      CACHE(2,NC) = XW1 - BNDRY(2)
      CACHE(3,NC) = BNDRY(3) - YW2
      CACHE(4,NC) = BNDRY(4) - XW2


      RETURN
      END


*=======================================================================
* MJD to/from Gregorian calendar date.
*
* Given:
*   CODE      I         If 1, compute MJD from IY,IM,ID, else vice
*                       versa.
*
* Given and/or returned:
*   MJD       D         Modified Julian date, (MJD = JD - 2400000.5).
*   IY        I         Year.
*   IM        I         Month (for CODE.EQ.1 may exceed 12, or be zero,
*                       or negative).
*   ID        I         Day of month (for CODE.EQ.1, may exceed the
*                       legitimate number of days in the month, or be
*                       zero, or negative).
*
* Notes:
*   These algorithms are from D.A. Hatcher, QJRAS 25, 53-55, as modified
*   by P.T. Wallace for use in SLALIB (subroutines CLDJ and DJCL).
*
*   Author: Mark Calabretta, Australia Telescope National Facility
*-----------------------------------------------------------------------
      SUBROUTINE PGMJD (CODE, MJD, IY, IM, ID)
*-----------------------------------------------------------------------
      INTEGER   CODE, DD, ID, IM, IY, JD, JM, JY, N4
      DOUBLE PRECISION MJD
*-----------------------------------------------------------------------
      IF (CODE.EQ.1) THEN
*       Calendar date to MJD.
        IF (IM.LT.1) THEN
          JY = IY - 1 + IM/12
          JM = 12 + MOD(IM,12)
        ELSE
          JY = IY + (IM-1)/12
          JM = MOD(IM-1,12) + 1
        END IF

        MJD =DBLE((1461*(JY - (12-JM)/10 + 4712))/4
     :          +(306*MOD(JM+9, 12) + 5)/10
     :          -(3*((JY - (12-JM)/10 + 4900)/100))/4
     :          + ID - 2399904)

      ELSE
*       MJD to calendar date.
        JD = 2400001 + INT(MJD)

        N4 =  4*(JD + ((2*((4*JD - 17918)/146097)*3)/4 + 1)/2 - 37)
        DD = 10*(MOD(N4-237, 1461)/4) + 5

        IY = N4/1461 - 4712
        IM = MOD(2 + DD/306, 12) + 1
        ID = MOD(DD, 306)/10 + 1
      END IF

      RETURN
      END



*=======================================================================
* This FORTRAN wrapper on PGSBOX exists solely to define fixed-length
* CHARACTER arguments for cpgsbox().
*
* Author: Mark Calabretta, Australia Telescope National Facility
*-----------------------------------------------------------------------
      SUBROUTINE PGSBOK (BLC, TRC, IDENTS, OPT, LABCTL, LABDEN, CI,
     :  GCODE, TIKLEN, NG1, GRID1, NG2, GRID2, DOEQ, NLFUNC, NLC, NLI,
     :  NLD, NLCPRM, NLIPRM, NLDPRM, NC, IC, CACHE, IERR)
*-----------------------------------------------------------------------
      LOGICAL   DOEQ
      INTEGER   CI(7), GCODE(2), IC, IERR, LABCTL, LABDEN, NC, NG1, NG2,
     :          NLC, NLD, NLI, NLIPRM(NLI)
      REAL      BLC(2), TRC(2)
      DOUBLE PRECISION CACHE(4,0:NC), GRID1(0:NG1), GRID2(0:NG2),
     :          NLDPRM(NLD), TIKLEN
      CHARACTER IDENTS(3)*80, NLCPRM(NLC)*1, OPT(2)*1

      EXTERNAL NLFUNC
*-----------------------------------------------------------------------
      CALL PGSBOX (BLC, TRC, IDENTS, OPT, LABCTL, LABDEN, CI, GCODE,
     :  TIKLEN, NG1, GRID1, NG2, GRID2, DOEQ, NLFUNC, NLC, NLI, NLD,
     :  NLCPRM, NLIPRM, NLDPRM, NC, IC, CACHE, IERR)

      RETURN
      END



*=======================================================================
* This FORTRAN wrapper on PGLBOX exists solely to define fixed-length
* CHARACTER arguments for cpglbox().
*
* Author: Mark Calabretta, Australia Telescope National Facility
*-----------------------------------------------------------------------
      SUBROUTINE PGLBOK (IDENTS, OPT, LABCTL, LABDEN, CI, GCODE, TIKLEN,
     :  NG1, GRID1, NG2, GRID2, DOEQ, NC, IC, CACHE, IERR)
*-----------------------------------------------------------------------
      LOGICAL   DOEQ
      INTEGER   CI(7), GCODE(2), IC, IERR, LABCTL, LABDEN, NC, NG1, NG2
      DOUBLE PRECISION CACHE(4,0:NC), GRID1(0:NG1), GRID2(0:NG2), TIKLEN
      CHARACTER IDENTS(3)*80, OPT(2)*1
*-----------------------------------------------------------------------
      CALL PGLBOX (IDENTS, OPT, LABCTL, LABDEN, CI, GCODE, TIKLEN, NG1,
     :  GRID1, NG2, GRID2, DOEQ, NC, IC, CACHE, IERR)

      RETURN
      END