C
C $Id: EzmapDemo.f,v 1.10 2005/06/22 21:44:24 kennison Exp $
C                                                                      
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C This file is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published
C by the Free Software Foundation; either version 2 of the License, or
C (at your option) any later version.
C
C This software is distributed in the hope that it will be useful, but
C WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
C General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this software; if not, write to the Free Software
C Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
C USA.
C

      PROGRAM EMDEMO
C
C Declarations required to do color fill.
C
        PARAMETER (LAMA=2000000,NCRA=100000,NGPS=10,LRWK=2*NCRA)
C
        DIMENSION IAMA(LAMA),XCRA(NCRA),YCRA(NCRA),IAAI(NGPS),IAGI(NGPS)
C
        DIMENSION RWRK(LRWK)
        EQUIVALENCE (RWRK(1),XCRA(1)),(RWRK(NCRA+1),YCRA(1))
C
        EXTERNAL COLORA,DRAWLA,COLORB,DRAWLB
C
C Declare some internal common blocks so as to be able to retrieve
C various values for the State Plane projection.
C
        COMMON /MAPCM2/  BLAM,BLOM,PEPS,SLAM,SLOM,UCEN,UMAX,UMIN,UOFF,
     +                   URNG,VCEN,VMAX,VMIN,VOFF,VRNG,ISSL
        DOUBLE PRECISION BLAM,BLOM,PEPS,SLAM,SLOM,UCEN,UMAX,UMIN,UOFF,
     +                   URNG,VCEN,VMAX,VMIN,VOFF,VRNG
        INTEGER          ISSL
        SAVE   /MAPCM2/
C
        COMMON /USGSC1/ UTPA(15),UUMN,UUMX,UVMN,UVMX,IPRF
          DOUBLE PRECISION UTPA,UUMN,UUMX,UVMN,UVMX
          INTEGER IPRF
        SAVE   /USGSC1/
C
C Declare a common block in which the directory name will be passed.
C
        COMMON /EDRGDI/ USNM
          CHARACTER*60 USNM
        SAVE   /EDRGDI/
C
C Declare a couple of variables that determine what projection will be
C used and how the portion of the projection to be shown is selected.
C
        CHARACTER*2 PTYP,LTYP,CTMP
C
C Declare a variable to determine what outline dataset will be used.
C
        CHARACTER*2 ODNM
C
C Declare some two-element arrays to hold values determining what
C portion of the projection is to be shown.
C
        DIMENSION PLM1(2),PLM2(2),PLM3(2),PLM4(2)
C
C Declare a character variable into which to read commands.
C
        CHARACTER*1 COMD
C
C Declare a character variable to hold a plot label.
C
        CHARACTER*60 LABL
C
C Declare a temporary character variable to read labels and names into.
C
        CHARACTER*60 CHRT
C
C Declare the function and variable necessary to deal with area names.
C
        CHARACTER*128 MDFNME,FNME
C
C Declare a couple of arrays in which to define boxes to be filled.
C
        DIMENSION XBOX(4),YBOX(4)
C
C Declare the double-precision parameter array required for initializing
C the USGS transformations.
C
        DOUBLE PRECISION PARA(15)
C
C Declare arrays in which to put State Plane zone-number information.
C
        DIMENSION IZN0(139),IZN8(139)
C
C Declare arrays in which to save the coordinates of zone boundaries.
C
        DIMENSION PLT0(200,139),PLN0(200,139)
        DIMENSION PLT8(200,139),PLN8(200,139)
C
C Declare arrays in which to save the coordinates of zone center points
C (used to position the labels).
C
        DIMENSION CLT0(139),CLN0(139),CLT8(139),CLN8(139)
C
C Declare some double precision temporaries to be used in generating
C the Piet Hein super-ellipses that are used for the zone borders in
C place of rectangles.
C
        DOUBLE PRECISION DBLP,DANG,DSNA,DCSA
C
C Declare arrays to hold color indices to be passed to MDRGSC, defining
C the colors to be used on maps drawn from the RANGS/GSHHS database.
C
        DIMENSION LCOL(5),LCSF(5)
C
C Set the default value of the parameter that says whether or not the
C X workstation will be updated as changes are made.  (When making a
C movie, it's a good idea to set this to zero.)
C
        DATA IUXW / 1 /
C
C Set the default value of the flag that says whether the basic color
C scheme will be white on black (IBOW=0) or black on white (IBOW=1).
C
        DATA IBOW / 0 /
C
C Declare the "style" of plot to be done.
C
        DATA ISTY / 1 /
C
C Set the default values of the variables that position the map on the
C plotter frame.
C
        DATA XVPL,XVPR,YVPB,YVPT / .1,.9,.1,.9 /
C
C Set the default values of the parameters that determine the projection
C to be used.
C
        DATA PTYP,PLAT,PLON,ROTA / 'CE',3*0. /
C
C Set the default values of the parameters that determine what portion
C of the projection is shown.
C
        DATA LTYP,PLM1,PLM2,PLM3,PLM4 / 'MA',8*0. /
C
C Set the default value of the satellite altitude and the angles that
C determine where its camera is pointed.
C
        DATA SALT,SAN1,SAN2 / 3*0. /
C
C Set the default value of the parameter that determines what outline
C dataset will be used.
C
        DATA ODNM / 'CO' /
C
C Set the default value of the parameter saying at what level the new
C database "Earth..1" should be used.
C
        DATA ILVL / 3 /
C
C Set the default value of the parameter saying at what level the new
C RANGS/GSHHS database should be used.
C
        DATA IRGL / 5 /
C
C Set the default value of the miscellaneous minor parameters.
C
        DATA IPER,IELL,ILBL,IDOT,IDBD,GSPC,GSLA,GSLO,GLAT,GLON
     +     /    1,   0,   1,   0,  96,  5.,  0.,  0.,  0., 90. /
C
        DATA IINT,IRPF
     +     /   50,   1 /
C
C Set a couple of counters that keep track of how many NCGM and
C PostScript frames have been saved.
C
        DATA NFNG,NFPS / 0 , 0 /
C
C Set the default values of the parameters required to initialize the
C USGS transformations.
C
        DATA IPRJ,IZON,ISPH,PARA / 2 , 4400 , 8 , 15*0.D0 /
C
C Set the default value of the index that can be used to select a
C pre-defined plot.
C
        DATA IPDP / 1 /
C
C Set the default value of the flag that says whether to draw an
C informational label at the bottom of the plot.
C
        DATA ILPT / 1 /
C
C Set the default value of the major label at the top of the plot.
C
        DATA LABL / ' ' /
C
C Define the flag that controls whether or not lat/lon labels are to
C be displayed and the variable that specifies the width of characters
C to use in them.
C
        DATA ILLB,CSLL / -1,.012 /
C
C Define some quantities for the test during which the projection of
C a global spiral is drawn.
C
        DATA ISPC,SPLT,SPLN,NSPT / 5 , 30. , 30. , 400 /
C
C Define some quantities affecting the test during which fields of
C latitude and longitudes generated by using the inverse routines
C are contoured.
C
        DATA ICNC,NCDM / 6 , 501 /
C
C Define a constant that says how nearly rectangular the Piet Hein
C super-ellipses are to be.  (They are used in place of rectangles
C for the State Plane zones.)
C
        DATA DBLP / 10.D0 /
C
C Define a flag that says which set of State Plane zones to display.
C The value "-1" means "neither"; "0" means "NAD1927"; and "8" means
C "NAD1983".  Also define a variable that specifies the width of the
C characters to use in labeling the zones.
C
        DATA ISPG,CHSZ / -1 , .1 /
C
C Define the State Plane zone numbers.
C
        DATA (IZN0(I),IZN8(I),I=1,50)/
     +    0101,0101,0102,0102,5001,5001,5002,5002,5003,5003,
     +    5004,5004,5005,5005,5006,5006,5007,5007,5008,5008,
     +    5009,5009,5010,5010,0201,0201,0202,0202,0203,0203,
     +    0301,0301,0302,0302,0401,0401,0402,0402,0403,0403,
     +    0404,0404,0405,0405,0406,0406,0407,   0,0501,0501,
     +    0502,0502,0503,0503,0600,0600,0700,0700,1901,1901,
     +    0901,0901,0902,0902,0903,0903,1001,1001,1002,1002,
     +    5101,5101,5102,5102,5103,5103,5104,5104,5105,5105,
     +    1101,1101,1102,1102,1103,1103,1201,1201,1202,1202,
     +    1301,1301,1302,1302,1401,1401,1402,1402,1501,1501/
C
        DATA (IZN0(I),IZN8(I),I=51,100)/
     +    1502,1502,1601,1601,1602,1602,1701,1701,1702,1702,
     +    1703,1703,1801,1801,1802,1802,1900,1900,2001,2001,
     +    2002,2002,2101,   0,2102,   0,2103,   0,2111,2111,
     +    2112,2112,2113,2113,2201,2201,2202,2202,2203,2203,
     +    2301,2301,2302,2302,2401,2401,2402,2402,2403,2403,
     +       0,2500,2501,   0,2502,   0,2503,   0,   0,2600,
     +    2601,   0,2602,   0,2701,2701,2702,2702,2703,2703,
     +    2800,2800,2900,2900,3001,3001,3002,3002,3003,3003,
     +    3101,3101,3102,3102,3103,3103,3104,3104,3200,3200,
     +    3301,3301,3302,3302,3401,3401,3402,3402,3501,3501/
C
        DATA (IZN0(I),IZN8(I),I=101,139)/
     +    3502,3502,3601,3601,3602,3602,3701,3701,3702,3702,
     +    3800,3800,   0,3900,3901,   0,3902,   0,4001,4001,
     +    4002,4002,4100,4100,4201,4201,4202,4202,4203,4203,
     +    4204,4204,4205,4205,4301,4301,4302,4302,4303,4303,
     +    4400,4400,4501,4501,4502,4502,4601,4601,4602,4602,
     +    4701,4701,4702,4702,4801,4801,4802,4802,4803,4803,
     +    4901,4901,4902,4902,4903,4903,4904,4904,5201,5200,
     +    5201,5200,5202,5200,5300,   0,5400,   0/
C
C Open GKS.
C
        CALL GOPKS (6,0)
C
C Open and activate an X workstation.
C
        CALL GOPWK (1,0,8)
        CALL GACWK (1)
C
C Open and activate an NCGM workstation.
C
        CALL GOPWK (2,0,1)
        CALL GACWK (2)
C
C Open and activate a PostScript workstation.
C
        CALL GOPWK (3,0,20)
        CALL GACWK (3)
C
C Turn clipping off.
C
        CALL GSCLIP (0)
C
C Select a "fill area interior style" of "solid".
C
        CALL GSFAIS (1)
C
C Define colors on all the workstations.
C
        DO 101 IWID=1,3
          IF (IBOW.EQ.0) THEN
            CALL GSCR (IWID,  0,0.,0.,0.)
            CALL GSCR (IWID,  1,1.,1.,1.)
          ELSE
            CALL GSCR (IWID,  0,1.,1.,1.)
            CALL GSCR (IWID,  1,0.,0.,0.)
          END IF
          CALL GSCR (IWID,  2,1.,0.,0.)
          CALL GSCR (IWID,  3,0.,1.,0.)
          CALL GSCR (IWID,  4,0.,0.,1.)
          CALL GSCR (IWID,  5,0.,1.,1.)
          CALL GSCR (IWID,  6,1.,0.,1.)
          CALL GSCR (IWID,  7,1.,1.,0.)
          CALL GSCR (IWID,  8,.7,.7,0.)
          CALL GSCR (IWID,  9,1.,1.,0.)
          CALL GSCR (IWID, 10,0.,0.,0.)
          CALL GSCR (IWID, 11,1.,1.,.6)
          CALL GSCR (IWID, 12,1.,0.,0.)
          CALL GSCR (IWID,101,.2,.2,.8)
          CALL GSCR (IWID,102,.6,.4,.2)
          CALL GSCR (IWID,103,.2,.6,.4)
          CALL GSCR (IWID,104,.2,.8,.2)
          CALL GSCR (IWID,105,.4,.6,.2)
          CALL GSCR (IWID,106,.2,.4,.6)
          CALL GSCR (IWID,107,.6,.6,.6)
          CALL GSCR (IWID,201,0.,0.,1.)  !  RANGS/GSHHS - Ocean
          CALL GSCR (IWID,202,0.,1.,0.)  !  RANGS/GSHHS - Land
          CALL GSCR (IWID,203,0.,1.,1.)  !  RANGS/GSHHS - Lake
          CALL GSCR (IWID,204,0.,1.,0.)  !  RANGS/GSHHS - Island in lake
          CALL GSCR (IWID,205,0.,1.,1.)  !  RANGS/GSHHS - Pond on island
          CALL GSCR (IWID,206,1.,1.,1.)  !  RANGS/GSHHS - Outlines
  101   CONTINUE
C
C Deactivate the NCGM and PostScript workstations.
C
        CALL GDAWK (2)
        CALL GDAWK (3)
C
C Put SETER into recovery mode.  No real error recovery is done in
C this program, but attempts to draw a plot that lead to an error
C are detected and the user is given another chance to avoid the
C error by changing parameter values.
C
        CALL ENTSR (IELD,1)
C
C Initialize PLOTCHAR font, character color, outline color, and box
C fill parameters.
C
        CALL PCSETI ('FN - FONT NUMBER',25)
        CALL PCSETI ('CC - CHARACTER COLOR',11)
        CALL PCSETI ('OF - OUTLINE FLAG',1)
        CALL PCSETI ('OC - OUTLINE LINE COLOR',11)
        CALL PCSETI ('BF - BOX FLAG',3)
        CALL PCSETI ('BC(1) - BOX COLOR 1',12)
        CALL PCSETI ('BC(2) - BOX COLOR 2',10)
        CALL PCSETR ('BL - BOX LINE WIDTH',2.)
        CALL PCSETR ('BM - BOX MARGIN WIDTH',.5)
C
C Tell EZMAP what colors to use for the outlines and filled areas
C produced from the RANGS/GSHHS data.
C
        LCOL(1)=206
        LCOL(2)=206
        LCOL(3)=206
        LCOL(4)=206
        LCOL(5)=206
C
        LCSF(1)=201
        LCSF(2)=202
        LCSF(3)=203
        LCSF(4)=204
        LCSF(5)=205
C
        CALL MDRGSC (LCOL,LCSF)
C
C Generate boundary-line data for each of the State Plane zones.
C
        PRINT * , ' '
        PRINT * , 'Generating boundary-line data for State Plane zones.'
C
        CALL MAPROJ ('UT',0.,0.,0.)
C
        DO 104 ISPI=1,139
C
          IF (IZN0(ISPI).NE.0) THEN
C
            CALL ZERODA (PARA,15)
            CALL MPUTIN (2,IZN0(ISPI),0,PARA,0.D0,0.D0,0.D0,0.D0)
C
            EPSI=.001*(UMAX-UMIN)
C
            UMID=(UMIN+UMAX)/2.
            VMID=(VMIN+VMAX)/2.
            UDIF=.999*(UMAX-UMID)
            VDIF=.999*(VMAX-VMID)
C
            CALL MAPTRI (UMID,VMID,CLT0(ISPI),CLN0(ISPI))
C
            DO 102 I=1,200
C
              DANG=1.8D0*DBLE(I-1)*.017453292519943D0
              DSNA=SIN(DANG)
              DCSA=COS(DANG)
              R=1./REAL(((ABS(DCSA)/DBLE(UDIF))**DBLP+
     +                   (ABS(DSNA)/DBLE(VDIF))**DBLP)**(1.D0/DBLP))
              CALL MAPTRI (UMID+R*REAL(DCSA),VMID+R*REAL(DSNA),
     +                     PLT0(I,ISPI),PLN0(I,ISPI))
  102       CONTINUE
C
          END IF
C
          IF (IZN8(ISPI).NE.0) THEN
C
            CALL ZERODA (PARA,15)
            CALL MPUTIN (2,IZN8(ISPI),8,PARA,0.D0,0.D0,0.D0,0.D0)
C
            EPSI=.001*(UMAX-UMIN)
C
            UMID=(UMIN+UMAX)/2.
            VMID=(VMIN+VMAX)/2.
            UDIF=.999*(UMAX-UMID)
            VDIF=.999*(VMAX-VMID)
C
            CALL MAPTRI (UMID,VMID,CLT8(ISPI),CLN8(ISPI))
C
            DO 103 I=1,200
              DANG=1.8D0*DBLE(I-1)*.017453292519943D0
              DSNA=SIN(DANG)
              DCSA=COS(DANG)
              R=1./REAL(((ABS(DCSA)/DBLE(UDIF))**DBLP+
     +                   (ABS(DSNA)/DBLE(VDIF))**DBLP)**(1.D0/DBLP))
              CALL MAPTRI (UMID+R*REAL(DCSA),VMID+R*REAL(DSNA),
     +                     PLT8(I,ISPI),PLN8(I,ISPI))
  103       CONTINUE
C
          END IF
C
  104   CONTINUE
C
C Initialize the variable that will be returned to EZMAP when it asks
C for the name of the directory containing the RANGS/GSHHS data.
C
        USNM=' '
C
C If the X workstation is not to be updated, jump to get user input.
C
  105   IF (IUXW.EQ.0) GO TO 106
C
C Clear the X workstation.
C
        PRINT * , ' '
        PRINT * , 'Clearing the workstation.'
C
        CALL GCLRWK (1,1)
C
C Draw a map.
C
        PRINT * , ' '
        PRINT * , 'Drawing the map.'
C
        CALL MAPPOS (XVPL,XVPR,YVPB,YVPT)
        CALL MAPROJ (PTYP,PLAT,PLON,ROTA)
        IF (PTYP.EQ.'SV') THEN
          CALL MPSETR ('SA',SALT)
          CALL MPSETR ('S1',SAN1)
          CALL MPSETR ('S2',SAN2)
        END IF
        CALL MAPSET (LTYP,PLM1,PLM2,PLM3,PLM4)
        CALL MPSETI ('PE',IPER)
        CALL MPSETI ('EL',IELL)
        CALL MPSETI ('LA',ILBL)
        CALL MPSETI ('DO',IDOT)
        CALL MPSETI ('DD',IDBD)
        CALL MPSETR ('GR',GSPC)
        CALL MPSETR ('GT',GSLA)
        CALL MPSETR ('GN',GSLO)
        CALL MPSETR ('GP',1000.*GLAT+GLON)
        CALL MPSETI ('II',1000*IINT+IINT)
        CALL MPSETI ('RP',IRPF)
        IF (PTYP.NE.'UT') THEN
          CALL MAPINT
        ELSE
          CALL MPUTIN (IPRJ,IZON,ISPH,PARA,0.D0,0.D0,0.D0,0.D0)
        END IF
        IF (NERRO(NERR).NE.0) GO TO 901
        CALL GETSET (RVPL,RVPR,RVPB,RVPT,XWDL,XWDR,YWDB,YWDT,LNLG)
        EPSI=1.E-6*MIN(RVPR-RVPL,RVPT-RVPB)
        RVPL=RVPL+EPSI
        RVPR=RVPR-EPSI
        RVPB=RVPB+EPSI
        RVPT=RVPT-EPSI
        PRINT * , ' '
        PRINT * , 'U/V limits used: ', XWDL,XWDR,YWDB,YWDT
        IF (ODNM.EQ.'RG'.AND.IRGL.GE.5) THEN
          CALL MDRGDL (IRGL)
          IRGL=IRGL+5
          PRINT * , ' '
          PRINT * , 'RANGS/GSHHS data level picked:',MOD(IRGL,5)
        END IF
        IF (ISTY.EQ.0) THEN
          CALL MAPGRD
          CALL MAPLBL
          IF (ODNM.EQ.'E1'.OR.ODNM.EQ.'E2'.OR.ODNM.EQ.'E3') THEN
            IF (ODNM.EQ.'E1') CALL MPLNDR ('Earth..1',ILVL)
            IF (ODNM.EQ.'E2') CALL MPLNDR ('Earth..2',ILVL)
            IF (ODNM.EQ.'E3') CALL MPLNDR ('Earth..3',ILVL)
          ELSE IF (ODNM.EQ.'RG') THEN
            CALL MDRGOL (MOD(IRGL,5),RWRK,LRWK)
            IF (NERRO(NERR).NE.0) GO TO 901
          ELSE
            CALL MAPLOT
            IF (NERRO(NERR).NE.0) GO TO 901
          END IF
          CALL DRSPGD (ISPG,IZN0,PLT0,PLN0,CLT0,CLN0,
     +                      IZN8,PLT8,PLN8,CLT8,CLN8,CHSZ)
        ELSE
          IF (ODNM.EQ.'E1'.OR.ODNM.EQ.'E2'.OR.ODNM.EQ.'E3') THEN
            CALL MPSETI ('VS',20)
            CALL ARINAM (IAMA,LAMA)
            IF (NERRO(NERR).NE.0) GO TO 901
            IF (ODNM.EQ.'E1') CALL MPLNAM ('Earth..1',ILVL,IAMA)
            IF (ODNM.EQ.'E2') CALL MPLNAM ('Earth..2',ILVL,IAMA)
            IF (ODNM.EQ.'E3') CALL MPLNAM ('Earth..3',ILVL,IAMA)
            IF (NERRO(NERR).NE.0) GO TO 901
            CALL ARSCAM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,COLORB)
            IF (NERRO(NERR).NE.0) GO TO 901
            IF (ODNM.EQ.'E1') CALL MPLNDR ('Earth..1',ILVL)
            IF (ODNM.EQ.'E2') CALL MPLNDR ('Earth..2',ILVL)
            IF (ODNM.EQ.'E3') CALL MPLNDR ('Earth..3',ILVL)
            CALL MAPGRM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,DRAWLB)
            IF (NERRO(NERR).NE.0) GO TO 901
            IAIC=IGETAI (IAMA,CFUX(.5),CFUY(.5),IAAI,IAGI,NGPS,NOGU,1)
            PRINT * , ' '
            PRINT * , 'AREA IDENTIFIER AT CENTER: ',IAIC
            FNME=MDFNME(IAIC,1)
            PRINT * , 'NAME OF AREA AT CENTER: ',
     +                                   FNME(MDIFNB(FNME):MDILNB(FNME))
          ELSE IF (ODNM.EQ.'RG') THEN
            CALL MDRGSF (MOD(IRGL,5),RWRK,LRWK,IAMA,LAMA)
            IF (NERRO(NERR).NE.0) GO TO 901
            CALL MAPGRD
            IF (NERRO(NERR).NE.0) GO TO 901
          ELSE
            CALL MPSETC ('OU',ODNM)
            CALL MPSETI ('VS',0)
            CALL ARINAM (IAMA,LAMA)
            CALL MAPBLA (IAMA)
            IF (NERRO(NERR).NE.0) GO TO 901
            CALL ARSCAM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,COLORA)
            IF (NERRO(NERR).NE.0) GO TO 901
            CALL MAPLOT
            IF (NERRO(NERR).NE.0) GO TO 901
            CALL MAPGRM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,DRAWLA)
            IF (NERRO(NERR).NE.0) GO TO 901
            IAIC=IGETAI (IAMA,CFUX(.5),CFUY(.5),IAAI,IAGI,NGPS,NOGU,1)
            PRINT * , ' '
            PRINT * , 'AREA IDENTIFIER AT CENTER: ',IAIC
          END IF
          CALL MAPLMB
          CALL MAPLBL
          CALL DRSPGD (ISPG,IZN0,PLT0,PLN0,CLT0,CLN0,
     +                      IZN8,PLT8,PLN8,CLT8,CLN8,CHSZ)
        END IF
        IF (ILLB.NE.0) THEN
          P=REAL(SIGN(1,ILLB))
          CALL PCSETI ('BF - BOX FLAG',0)
          CALL MDLBLT (RVPL,RVPB,RVPL,RVPT,+P*CSLL,    0.,CSLL,0.,-P)
          CALL MDLBLT (RVPR,RVPB,RVPR,RVPT,-P*CSLL,    0.,CSLL,0.,+P)
          CALL MDLBLN (RVPL,RVPB,RVPR,RVPB,0.,+1.5*P*CSLL,CSLL,0.,0.)
          CALL MDLBLN (RVPL,RVPT,RVPR,RVPT,0.,-1.5*P*CSLL,CSLL,0.,0.)
          CALL PCSETI ('BF - BOX FLAG',3)
        END IF
        IF (LABL.NE.' ') CALL PLCHHQ (CFUX(.5),CFUY(.965),
     +                                LABL(1:MPILNB(LABL)),
     +                                          .018,0.,0.)
        IF (ILPT.NE.0) CALL DRPTLB (PTYP,IPRJ,IZON,ISPH,PARA,IPRF,
     +                                                  ODNM,ILVL)
C
C Flush buffers and update the X workstation.
C
        CALL SFLUSH
        CALL GUWK (1,0)
C
C Let the user enter a command.
C
  106   PRINT * , ' '
        PRINT * , 'Enter a command (H for help, Q for quit):'
C
        READ  '(A1)', COMD
C
        IF      (COMD.EQ.'H'.OR.COMD.EQ.'h') THEN
C
          PRINT * , ' '
          PRINT * , 'Possible commands are as follows:'
          PRINT * , ' '
          PRINT * , '  H  =>  Help (this is it)'
          PRINT * , '  C  =>  Change parameters (in several groups)'
          PRINT * , '  P  =>  Plot (specified predefined plots)'
          PRINT * , '  Q  =>  Quit (stop, halt, cease and desist)'
          PRINT * , '  R  =>  Redraw (for use after resizing X window)'
          PRINT * , '  S  =>  Save (to an NCGM or Postscript file)'
          PRINT * , '  T  =>  Tests (of three different kinds)'
          PRINT * , ' '
          PRINT * , 'In user responses, case is ignored.  When a yes'
          PRINT * , 'or no answer is requested, hitting RETURN means'
          PRINT * , 'no.  When a numeric value is requested, hitting'
          PRINT * , 'RETURN implies no change in the current value.'
          PRINT * , ' '
          PRINT * , 'Most EZMAP errors will not result in termination'
          PRINT * , 'of the demo.  It should usually be obvious what'
          PRINT * , 'needs to be changed to avoid the error.'
C
          GO TO 106
C
        ELSE IF (COMD.EQ.'C'.OR.COMD.EQ.'c') THEN
C
          PRINT * , ' '
          PRINT * , 'Changing internal parameter values.'
C
          PRINT * , ' '
          PRINT * , 'Change position of map on plotter frame (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
            PRINT * , 'Current left-edge value:  ',XVPL
            PRINT * , 'Current right-edge value: ',XVPR
            PRINT * , 'Current bottom-edge value:',YVPB
            PRINT * , 'Current top-edge value:   ',YVPT
            PRINT * , ' '
            PRINT * , 'Enter new left-edge value (between 0 and 1):'
            CALL EMRDRN (XVPL,XVPL)
            XVPL=MAX(0.,MIN(1.,XVPL))
            PRINT * , ' '
            PRINT * , 'Enter new right-edge value (between 0 and 1):'
            CALL EMRDRN (XVPR,XVPR)
            XVPR=MAX(0.,MIN(1.,XVPR))
            PRINT * , ' '
            PRINT * , 'Enter new bottom-edge value (between 0 and 1):'
            CALL EMRDRN (YVPB,YVPB)
            YVPB=MAX(0.,MIN(1.,YVPB))
            PRINT * , ' '
            PRINT * , 'Enter new top-edge value (between 0 and 1):'
            CALL EMRDRN (YVPT,YVPT)
            YVPT=MAX(0.,MIN(1.,YVPT))
          END IF
C
          PRINT * , ' '
          PRINT * , 'Change projection type, center, rotation angle,'
          PRINT * , 'other values for USGS transformations (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
            PRINT * , 'These are the possible projection types:'
            PRINT * , ' '
            PRINT * , '  LC => Lambert Conformal Conic'
            PRINT * , '  ST => Stereographic'
            PRINT * , '  OR => Orthographic'
            PRINT * , '  LE => Lambert Equal-area'
            PRINT * , '  GN => Gnomonic'
            PRINT * , '  AE => Azimuthal Equidistant'
            PRINT * , '  CE => Cylindrical Equidistant'
            PRINT * , '  EA => Cylindrical Equal-Area'
            PRINT * , '  ME => Mercator'
            PRINT * , '  RM => Rotated Mercator'
            PRINT * , '  MO => Mollweide'
            PRINT * , '  RO => Robinson'
            PRINT * , '  SV => Satellite-View'
            PRINT * , '  UT => USGS Transformations'
            PRINT * , ' '
            PRINT * , 'Current projection type: ',PTYP
            PRINT * , ' '
            PRINT * , 'Enter new projection type:'
            READ  '(A2)', CTMP
            CALL MUPPER (CTMP)
            IF (CTMP.EQ.'LC'.OR.CTMP.EQ.'ST'.OR.CTMP.EQ.'OR'.OR.
     +          CTMP.EQ.'LE'.OR.CTMP.EQ.'GN'.OR.CTMP.EQ.'AE'.OR.
     +          CTMP.EQ.'CE'.OR.CTMP.EQ.'EA'.OR.CTMP.EQ.'ME'.OR.
     +          CTMP.EQ.'RM'.OR.CTMP.EQ.'MO'.OR.CTMP.EQ.'RO'.OR.
     +          CTMP.EQ.'SV'.OR.CTMP.EQ.'UT') PTYP=CTMP
            IF (PTYP.EQ.'LC') THEN
              PRINT * , ' '
              PRINT * , 'Current central meridian:      ',PLON
              PRINT * , 'Current 1st standard parallel: ',PLAT
              PRINT * , 'Current 2nd standard parallel: ',ROTA
              PRINT * , ' '
              PRINT * , 'Enter new central meridian (between -180 and 18
     +0):'
              CALL EMRDRN (PLON,PLON)
              PLON=MAX(-180.,MIN(180.,PLON))
              PRINT * , ' '
              PRINT * , 'Enter new 1st standard parallel (between -90 an
     +d 90):'
              CALL EMRDRN (PLAT,PLAT)
              PLAT=MAX(-90.,MIN(90.,PLAT))
              PRINT * , ' '
              PRINT * , 'Enter new 2nd standard parallel (between -90 an
     +d 90):'
              CALL EMRDRN (ROTA,ROTA)
              ROTA=MAX(-90.,MIN(90.,ROTA))
            ELSE IF (PTYP.EQ.'RM') THEN
              PRINT * , ' '
              PRINT * , 'Current central meridian:      ',PLON
              PRINT * , 'Current rotation angle:        ',ROTA
              PRINT * , ' '
              PRINT * , 'Enter new central meridian (between -180 and 18
     +0):'
              CALL EMRDRN (PLON,PLON)
              PLON=MAX(-180.,MIN(180.,PLON))
              PRINT * , ' '
              PRINT * , 'Enter new rotation angle (between -180 and 180)
     +:'
              CALL EMRDRN (ROTA,ROTA)
              ROTA=MAX(-180.,MIN(180.,ROTA))
            ELSE IF (PTYP.EQ.'UT') THEN
              PRINT * , ' '
              PRINT * , 'For most of the USGS transformations, one may '
              PRINT * , 'select any one of various "spheroids"; each of'
              PRINT * , 'these is a particular ellipsoid having semi-  '
              PRINT * , 'major and semiminor axes of specified lengths '
              PRINT * , '(in meters).  A couple of the transformations '
              PRINT * , 'require using a particular spheroid and some  '
              PRINT * , 'require a "sphere of reference" with a radius '
              PRINT * , 'equal to the semimajor axis of the selected   '
              PRINT * , 'spheroid.  The possible spheroids are these:  '
              PRINT * , '                                              '
              PRINT * , '  -1 => arbitrary, user-specified             '
              PRINT * , '   0 => Clarke 1866                           '
              PRINT * , '   1 => Clarke 1880                           '
              PRINT * , '   2 => Bessel                                '
              PRINT * , '   3 => New International 1967                '
              PRINT * , '   4 => International 1909                    '
              PRINT * , '   5 => WGS 72                                '
              PRINT * , '   6 => Everest                               '
              PRINT * , '   7 => WGS 66                                '
              PRINT * , '   8 => GRS 1980                              '
              PRINT * , '   9 => Airy                                  '
              PRINT * , '  10 => Modified Everest                      '
              PRINT * , '  11 => Modified Airy                         '
              PRINT * , '  12 => WGS 84                                '
              PRINT * , '  13 => Southeast Asia                        '
              PRINT * , '  14 => Australian National                   '
              PRINT * , '  15 => Krassovsky                            '
              PRINT * , '  16 => Hough                                 '
              PRINT * , '  17 => Mercury 1960                          '
              PRINT * , '  18 => Modified Mercury 1968                 '
              PRINT * , '  19 => Normal Sphere                         '
              PRINT * , ' '
              PRINT * , 'The spheroid currently selected:',ISPH
              PRINT * , ' '
              PRINT * , 'Enter the desired spheroid number:'
              CALL EMRDIN (ISPH,ISPH)
              ISPH=MAX(-1,MIN(19,ISPH))
              IF (ISPH.LT.0) THEN
                PRINT * , ' '
                PRINT * , 'We need the length of the semimajor axis of '
                PRINT * , 'the ellipsoid (the equatorial radius), in   '
                PRINT * , 'meters.  The value must be greater than or  '
                PRINT * , 'equal to 1. (The actual equatorial radius   '
                PRINT * , 'of the earth is about 6,378,000 meters, but '
                PRINT * , 'a smaller value can be used if the scale of '
                PRINT * , 'the map is unimportant.)                    '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(1)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(1)))
                PARA(1)=DBLE(MAX(.1,TEMP))
                PRINT * , ' '
                PRINT * , 'We need the value of the square of the      '
                PRINT * , 'eccentricy of the ellipsoid; it must be     '
                PRINT * , 'less than .1 (generally, it''s much less).  '
                PRINT * , 'A value greater than 1 is treated as the    '
                PRINT * , 'length of the semiminor axis (the polar     '
                PRINT * , 'radius), in meters.  (The actual polar      '
                PRINT * , 'radius of the earth is about 6,357,000      '
                PRINT * , 'meters.)                                    '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(2)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(2)))
                PARA(2)=DBLE(MAX(0.,TEMP))
              END IF
              PRINT * , ' '
              PRINT * , 'Possible USGS transformations are as follows: '
              PRINT * , '                                              '
              PRINT * , '   1 => Universal Transverse Mercator (UTM)   '
              PRINT * , '   2 => State Plane Coordinates               '
              PRINT * , '   3 => Alber''s Equal-Area Conic             '
              PRINT * , '   4 => Lambert Conformal Conic               '
              PRINT * , '   5 => Mercator                              '
              PRINT * , '   6 => Polar Stereographic                   '
              PRINT * , '   7 => Polyconic                             '
              PRINT * , '   8 => Equidistant Conic                     '
              PRINT * , '   9 => Transverse Mercator                   '
              PRINT * , '  10 => Stereographic                         '
              PRINT * , '  11 => Lambert Azimuthal Equal-Area          '
              PRINT * , '  12 => Azimuthal Equidistant                 '
              PRINT * , '  13 => Gnomonic                              '
              PRINT * , '  14 => Orthographic                          '
              PRINT * , '  15 => Perspective                           '
              PRINT * , '  16 => Sinusoidal                            '
              PRINT * , '  17 => Equirectangular                       '
              PRINT * , '  18 => Miller Cylindrical                    '
              PRINT * , '  19 => Van der Grinten I                     '
              PRINT * , '  20 => Oblique Mercator                      '
              PRINT * , '  21 => Robinson                              '
              PRINT * , '  22 => Space Oblique Mercator                '
              PRINT * , '  23 => Modified Stereographic for Alaska     '
              PRINT * , ' '
              PRINT * , 'The transformation currently selected:',IPRJ
              PRINT * , ' '
              PRINT * , 'Enter the desired transformation number:'
              CALL EMRDIN (IPRJ,IPRJ)
              IPRJ=MAX(1,MIN(23,IPRJ))
              IF      (IPRJ.EQ.1) THEN  !  UTM Coordinates
                PRINT * , ' '
                PRINT * , 'The UTM system requires a zone number with  '
                PRINT * , 'an absolute value from 1 to 60, inclusive.  '
                PRINT * , 'Each zone is a strip 6 degrees of longitude '
                PRINT * , 'in width.  Zone 1 is from 180W to 174W, zone'
                PRINT * , '2 from 174W to 168W, and so on; zone 60 is  '
                PRINT * , 'from 174E to 180E.  Positive zone numbers   '
                PRINT * , 'imply the Northern hemisphere and negative  '
                PRINT * , 'ones the Southern hemisphere.               '
                PRINT * , ' '
                PRINT * , 'The zone number currently selected:',IZON
                PRINT * , ' '
                PRINT * , 'Enter the desired zone number:'
                CALL EMRDIN (IZON,IZON)
                IZON=MAX(-60,MIN(60,IZON))
                IF (IZON.EQ.0) IZON=13
              ELSE IF (IPRJ.EQ.2) THEN  !  State Plane Coordinates
                PRINT * , ' '
                PRINT * , 'The older version of the State Plane system '
                PRINT * , 'is based on the North American Datum of 1927'
                PRINT * , '(NAD1927) and the newer one is based on the '
                PRINT * , 'North American Datum of 1983 (NAD1983).  The'
                PRINT * , 'first is selected by using a spheroid with  '
                PRINT * , 'index 0 and the second is selected by using '
                PRINT * , 'a spheroid with index 8.                    '
                PRINT * , ' '
                PRINT * , 'The current spheroid index:',ISPH
                PRINT * , ' '
                PRINT * , 'Enter the desired spheroid index:'
                CALL EMRDIN (ISPH,ISPH)
                IF (ISPH.NE.0.AND.ISPH.NE.8) THEN
                  PRINT * , ' '
                  PRINT * , 'The value entered is invalid; spheroid 8  '
                  PRINT * , 'will be used.                             '
                  ISPH=8
                END IF
                PRINT * , ' '
                PRINT * , 'Each version of the State Plane system uses '
                PRINT * , 'zone numbers to identify a particular state '
                PRINT * , 'or portion of a state.  The zone numbers are'
                PRINT * , 'four digits long; values may be found in the'
                PRINT * , 'on-line documentation and on some predefined'
                PRINT * , 'plots available from this demo program.     '
                PRINT * , ' '
                PRINT * , 'The current zone number:',IZON
                PRINT * , ' '
                PRINT * , 'Enter the desired zone number:'
                CALL EMRDIN (IZON,IZON)
                IF (ISPH.EQ.0) THEN
                  DO 114 I=1,139
                    IF (IZON.EQ.IZN0(I)) GO TO 116
  114             CONTINUE
                ELSE
                  DO 115 I=1,139
                    IF (IZON.EQ.IZN8(I)) GO TO 116
  115             CONTINUE
                END IF
                PRINT * , ' '
                PRINT * , 'The selected zone number is invalid.  The   '
                PRINT * , 'value 2800 (New Hampshire) will be used as  '
                PRINT * , 'an illustration.                            '
                IZON=2800
  116           CONTINUE
              ELSE IF (IPRJ.EQ.3) THEN  !  Alber's Equal-Area Conic
                PRINT * , ' '
                PRINT * , 'We need the latitude of the first standard  '
                PRINT * , 'parallel (between -90 and +90).             '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(3)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(3)))
                PARA(3)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the second standard '
                PRINT * , 'parallel (between -90 and +90).             '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(4)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(4)))
                PARA(4)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the origin of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.4) THEN  !  Lambert Conformal Conic
                PRINT * , ' '
                PRINT * , 'We need the latitude of the first standard  '
                PRINT * , 'parallel (between -90 and +90).             '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(3)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(3)))
                PARA(3)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the second standard '
                PRINT * , 'parallel (between -90 and +90).             '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(4)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(4)))
                PARA(4)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the origin of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.5) THEN  !  Mercator
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of true scale (between '
                PRINT * , '-90 and +90).                               '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.6) THEN  !  Polar Stereographic
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of true scale (between '
                PRINT * , '-90 and +90).                               '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.7) THEN  !  Polyconic
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the origin of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.8) THEN  !  Equidistant Conic
                PRINT * , ' '
                PRINT * , 'This transformation has a flag which is 0   '
                PRINT * , 'if there is only one standard parallel and  '
                PRINT * , '1 if there are two standard parallels.      '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',INT(PARA(9)
     +)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDIN (ITMP,INT(PARA(9)))
                PARA(9)=DBLE(MAX(0,MIN(1,ITMP)))
                PRINT * , ' '
                IF (PARA(9).EQ.0.D0) THEN
                  PRINT * , 'We need the latitude of the standard      '
                  PRINT * , 'parallel (between -90 and +90).           '
                ELSE
                  PRINT * , 'We need the latitude of the first         '
                  PRINT * , 'standard parallel (between -90 and +90).  '
                END IF
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(3)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(3)))
                PARA(3)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                IF (PARA(9).NE.0.D0) THEN
                  PRINT * , ' '
                  PRINT * , 'We need the latitude of the second        '
                  PRINT * , 'standard parallel (between -90 and +90).  '
                  PRINT * , ' '
                  PRINT * , 'The current value is as follows:',PARA(4)
                  PRINT * , ' '
                  PRINT * , 'Enter the desired new value:'
                  CALL EMRDRN (TEMP,REAL(PARA(4)))
                  PARA(4)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                END IF
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the origin of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.9) THEN  !  Transverse Mercator
                PRINT * , ' '
                PRINT * , 'We need the scale factor at the central     '
                PRINT * , 'meridian (use a value close to 1).          '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(3)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(3)))
                PARA(3)=DBLE(MAX(.1,MIN(10.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the origin (between '
                PRINT * , '-90 and +90).                               '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.10) THEN  !  Stereographic
                PRINT * , ' '
                PRINT * , 'The Stereographic projection uses a sphere  '
                PRINT * , 'of reference, rather than an ellipsoid.     '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the center of the  '
                PRINT * , 'projection (between -180 and +180).         '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the center of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.11) THEN  !  Lambert Azimuthal Equal-Area
                PRINT * , ' '
                PRINT * , 'The Lambert Azimuthal Equal-Area projection '
                PRINT * , 'uses a sphere of reference, rather than an  '
                PRINT * , 'ellipsoid.                                  '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the center of the  '
                PRINT * , 'projection (between -180 and +180).         '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the center of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.12) THEN  !  Azimuthal Equidistant
                PRINT * , ' '
                PRINT * , 'The Azimuthal Equidistant projection uses   '
                PRINT * , 'a sphere of reference, rather than an       '
                PRINT * , 'ellipsoid.                                  '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the center of the  '
                PRINT * , 'projection (between -180 and +180).         '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the center of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.13) THEN  !  Gnomonic
                PRINT * , ' '
                PRINT * , 'The Gnomonic projection uses a sphere of    '
                PRINT * , 'reference, rather than an ellipsoid.        '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the center of the  '
                PRINT * , 'projection (between -180 and +180).         '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the center of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.14) THEN  !  Orthographic
                PRINT * , ' '
                PRINT * , 'The Orthographic projection uses a sphere of'
                PRINT * , 'reference, rather than an ellipsoid.        '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the center of the  '
                PRINT * , 'projection (between -180 and +180).         '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the center of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.15) THEN  !  Perspective
                PRINT * , ' '
                PRINT * , 'The Perspective projection uses a sphere of '
                PRINT * , 'reference, rather than an ellipsoid.        '
                PRINT * , ' '
                PRINT * , 'We need the height of the perspective point '
                PRINT * , 'above the surface (in meters).              '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(3)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(3)))
                PARA(3)=DBLE(MAX(0.,TEMP))
                PRINT * , ' '
                PRINT * , 'We need the longitude of the center of the  '
                PRINT * , 'projection (between -180 and +180).         '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of the center of the   '
                PRINT * , 'projection (between -90 and +90).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.16) THEN  !  Sinusoidal
                PRINT * , ' '
                PRINT * , 'The Sinusoidal projection uses a sphere of  '
                PRINT * , 'reference, rather than an ellipsoid.        '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.17) THEN  !  Equirectangular
                PRINT * , ' '
                PRINT * , 'The Equirectangular projection uses a       '
                PRINT * , 'sphere of reference, rather than an         '
                PRINT * , 'ellipsoid.                                  '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need the latitude of true scale (between '
                PRINT * , '-90 and +90).                               '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(6)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(6)))
                PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.18) THEN  !  Miller Cylindrical
                PRINT * , ' '
                PRINT * , 'The Miller Cylindrical projection uses a    '
                PRINT * , 'sphere of reference, rather than an         '
                PRINT * , 'ellipsoid.                                  '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.19) THEN  !  Van der Grinten I
                PRINT * , ' '
                PRINT * , 'The Van der Grinten projection uses a sphere'
                PRINT * , 'of reference, rather than an ellipsoid.     '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.20) THEN  !  Oblique Mercator
                PRINT * , ' '
                PRINT * , 'The Oblique Mercator transformation has a   '
                PRINT * , 'flag which is 0 if the azimuth angle is to  '
                PRINT * , 'be specified indirectly and 1 if it is to   '
                PRINT * , 'be specified directly.                      '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',INT(PARA(13
     +))
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDIN (ITMP,INT(PARA(13)))
                PARA(13)=DBLE(MAX(0,MIN(1,ITMP)))
                PRINT * , ' '
                PRINT * , 'We need the scale factor at the center of   '
                PRINT * , 'the projection.                             '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(3)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(3)))
                PARA(3)=DBLE(MAX(.1,MIN(10.,TEMP)))
                IF (PARA(13).EQ.0.D0) THEN
                  PRINT * , ' '
                  PRINT * , 'We need the latitude of the origin of the '
                  PRINT * , 'projection (between -90 and +90).         '
                  PRINT * , ' '
                  PRINT * , 'The current value is as follows:',PARA(6)
                  PRINT * , ' '
                  PRINT * , 'Enter the desired new value:'
                  CALL EMRDRN (TEMP,REAL(PARA(6)))
                  PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                  PRINT * , ' '
                  PRINT * , 'We need the longitude of point 1 on a line'
                  PRINT * , 'defining the azimuth angle (between -180  '
                  PRINT * , 'and +180).                                '
                  PRINT * , ' '
                  PRINT * , 'The current value is as follows:',PARA( 9)
                  PRINT * , ' '
                  PRINT * , 'Enter the desired new value:'
                  CALL EMRDRN (TEMP,REAL(PARA( 9)))
                  PARA( 9)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                  PRINT * , ' '
                  PRINT * , 'We need the latitude of point 1 on a line '
                  PRINT * , 'defining the azimuth angle (between -90   '
                  PRINT * , 'and +90).                                 '
                  PRINT * , ' '
                  PRINT * , 'The current value is as follows:',PARA(10)
                  PRINT * , ' '
                  PRINT * , 'Enter the desired new value:'
                  CALL EMRDRN (TEMP,REAL(PARA(10)))
                  PARA(10)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                  PRINT * , ' '
                  PRINT * , 'We need the longitude of point 2 on a line'
                  PRINT * , 'defining the azimuth angle (between -180  '
                  PRINT * , 'and +180).                                '
                  PRINT * , ' '
                  PRINT * , 'The current value is as follows:',PARA(11)
                  PRINT * , ' '
                  PRINT * , 'Enter the desired new value:'
                  CALL EMRDRN (TEMP,REAL(PARA(11)))
                  PARA(11)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                  PRINT * , ' '
                  PRINT * , 'We need the latitude of point 2 on a line '
                  PRINT * , 'defining the azimuth angle (between -90   '
                  PRINT * , 'and +90).                                 '
                  PRINT * , ' '
                  PRINT * , 'The current value is as follows:',PARA(12)
                  PRINT * , ' '
                  PRINT * , 'Enter the desired new value:'
                  CALL EMRDRN (TEMP,REAL(PARA(12)))
                  PARA(12)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                ELSE
                  PRINT * , ' '
                  PRINT * , 'We need the azimuth angle of the center   '
                  PRINT * , 'line, measured east of north (between -180'
                  PRINT * , 'and +180).                                '
                  PRINT * , ' '
                  PRINT * , 'The current value is as follows:',PARA(4)
                  PRINT * , ' '
                  PRINT * , 'Enter the desired new value:'
                  CALL EMRDRN (TEMP,REAL(PARA(4)))
                  PARA(4)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                  PRINT * , ' '
                  PRINT * , 'We need the longitude of the origin of the'
                  PRINT * , 'projection (between -180 and +180).       '
                  PRINT * , ' '
                  PRINT * , 'The current value is as follows:',PARA(5)
                  PRINT * , ' '
                  PRINT * , 'Enter the desired new value:'
                  CALL EMRDRN (TEMP,REAL(PARA(5)))
                  PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                  PRINT * , ' '
                  PRINT * , 'We need the latitude of the origin of the '
                  PRINT * , 'projection (between -90 and +90).         '
                  PRINT * , ' '
                  PRINT * , 'The current value is as follows:',PARA(6)
                  PRINT * , ' '
                  PRINT * , 'Enter the desired new value:'
                  CALL EMRDRN (TEMP,REAL(PARA(6)))
                  PARA(6)=DBLE(MAX(-90.,MIN(90.,TEMP)))
                END IF
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.21) THEN  !  Robinson
                PRINT * , ' '
                PRINT * , 'The Robinson projection uses a sphere of    '
                PRINT * , 'reference, rather than an ellipsoid.  This  '
                PRINT * , 'projection is much used by the National     '
                PRINT * , 'Geographic Society for world maps.          '
                PRINT * , ' '
                PRINT * , 'We need the longitude of the central        '
                PRINT * , 'meridian (between -180 and +180).           '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(5)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(5)))
                PARA(5)=DBLE(MAX(-180.,MIN(180.,TEMP)))
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.22) THEN  !  Space Oblique Mercator
                PRINT * , ' '
                PRINT * , 'The Space Oblique Mercator appears to be    '
                PRINT * , 'set up so that one can specify a particular '
                PRINT * , 'Landsat series and a particular orbit and   '
                PRINT * , 'get a map that shows with minimal distortion'
                PRINT * , 'the part of the earth passed over during    '
                PRINT * , 'that particular orbit.                      '
                PRINT * , ' '
                PRINT * , 'We need a Landsat number (an integer between'
                PRINT * , '1 and 5, inclusive).                        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',INT(PARA(3)
     +)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDIN (ITMP,INT(PARA(3)))
                PARA(3)=DBLE(MAX(1,MIN(5,ITMP)))
                PRINT * , ' '
                PRINT * , 'We need a path number (an integer between 1 '
                PRINT * , 'and 251, inclusive, for Landsat numbers 1, 2'
                PRINT * , 'and 3, and between 1 and 233, inclusive, for'
                PRINT * , 'Landsat numbers 4 and 5).                   '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',INT(PARA(4)
     +)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDIN (ITMP,INT(PARA(4)))
                IF (INT(PARA(3)).LE.3) THEN
                  PARA(4)=DBLE(MAX(1,MIN(251,ITMP)))
                ELSE
                  PARA(4)=DBLE(MAX(1,MIN(233,ITMP)))
                END IF
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE IF (IPRJ.EQ.23) THEN  !  Modified Stereographic (AK)
                PRINT * , ' '
                PRINT * , 'The Modified Stereographic projection for   '
                PRINT * , 'Alaska uses the Clarke 1866 spheroid, no    '
                PRINT * , 'matter what.                                '
                PRINT * , ' '
                PRINT * , 'We need a false easting (in meters).        '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(7)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(7)))
                PARA(7)=DBLE(TEMP)
                PRINT * , ' '
                PRINT * , 'We need a false northing (in meters).       '
                PRINT * , ' '
                PRINT * , 'The current value is as follows:',PARA(8)
                PRINT * , ' '
                PRINT * , 'Enter the desired new value:'
                CALL EMRDRN (TEMP,REAL(PARA(8)))
                PARA(8)=DBLE(TEMP)
              ELSE
                IZON=0
              END IF
            ELSE
              PRINT * , ' '
              PRINT * , 'Current pole latitude:  ',PLAT
              PRINT * , 'Current pole longitude: ',PLON
              PRINT * , 'Current rotation angle: ',ROTA
              PRINT * , ' '
              PRINT * , 'Enter new pole latitude (between -90 and 90):'
              CALL EMRDRN (PLAT,PLAT)
              PLAT=MAX(-90.,MIN(90.,PLAT))
              PRINT * , ' '
              PRINT * , 'Enter new pole longitude (between -180 and 180)
     +:'
              CALL EMRDRN (PLON,PLON)
              PLON=MAX(-180.,MIN(180.,PLON))
              PRINT * , ' '
              PRINT * , 'Enter new rotation angle (between -180 and 180)
     +:'
              CALL EMRDRN (ROTA,ROTA)
              ROTA=MAX(-180.,MIN(180.,ROTA))
            END IF
            IF (PTYP.EQ.'SV') THEN
              PRINT * , ' '
              PRINT * , 'Satellite-view projection will be used.'
              PRINT * , ' '
              PRINT * , 'Current satellite height :',SALT
              PRINT * , 'Current satellite angle 1:',SAN1
              PRINT * , 'Current satellite angle 2:',SAN2
              PRINT * , ' '
              PRINT * , 'Enter new satellite height (measured from'
              PRINT * , 'earth center, 1 or greater):'
              CALL EMRDRN (SALT,SALT)
              IF (SALT.LE.1.) SALT=0.
              PRINT * , ' '
              PRINT * , 'Enter new satellite angle 1 (0 to 180):'
              CALL EMRDRN (SAN1,SAN1)
              SAN1=MAX(0.,MIN(180.,SAN1))
              PRINT * , ' '
              PRINT * , 'Enter new satellite angle 2 (-180 to 180):'
              CALL EMRDRN (SAN2,SAN2)
              SAN2=MAX(-180.,MIN(180.,SAN2))
            END IF
          END IF
C
          PRINT * , ' '
          PRINT * , 'Change limits selecting part of projection shown (Y
     + or N)?'
          READ  '(A1)', COMD
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
  108       PRINT * , ' '
            PRINT * , 'Possible limit selection types:'
            PRINT * , ' '
            PRINT * , '  MA => Maximum useful area'
            PRINT * , '  CO => Corners opposite one another'
            PRINT * , '  PO => Points on the edges'
            PRINT * , '  AN => Angular extents'
            PRINT * , '  LI => Limits in the U/V plane'
            PRINT * , '  GR => Grid of lat/lon values'
            PRINT * , ' '
            PRINT * , 'Current limit selection type: ',LTYP
            PRINT * , ' '
            PRINT * , 'Enter new limit selection type:'
            READ  '(A2)', CTMP
            CALL MUPPER (CTMP)
            IF (PTYP.EQ.'LC'.AND.CTMP.EQ.'AN') THEN
              PRINT * , ' '
              PRINT * , 'Can''t use angular extents with a Lambert confo
     +rmal conic projection.'
              GO TO 108
            END IF
            IF (CTMP.EQ.'MA'.OR.CTMP.EQ.'CO'.OR.CTMP.EQ.'PO'.OR.
     +          CTMP.EQ.'AN'.OR.CTMP.EQ.'LI'.OR.CTMP.EQ.'GR') LTYP=CTMP
            IF (LTYP.EQ.'CO') THEN
              PRINT * , ' '
              PRINT * , 'Current latitude, 1st corner point: ',PLM1(1)
              PRINT * , 'Current longitude, 1st corner point:',PLM2(1)
              PRINT * , 'Current latitude, 2nd corner point: ',PLM3(1)
              PRINT * , 'Current longitude, 2nd corner point:',PLM4(1)
              PRINT * , ' '
              PRINT * , 'Enter new latitude of 1st corner point:'
              CALL EMRDRN (PLM1(1),PLM1(1))
              PRINT * , ' '
              PRINT * , 'Enter new longitude of 1st corner point:'
              CALL EMRDRN (PLM2(1),PLM2(1))
              PRINT * , ' '
              PRINT * , 'Enter new latitude of 2nd corner point:'
              CALL EMRDRN (PLM3(1),PLM3(1))
              PRINT * , ' '
              PRINT * , 'Enter new longitude of 2nd corner point:'
              CALL EMRDRN (PLM4(1),PLM4(1))
              PLM1(2)=0.
              PLM2(2)=0.
              PLM3(2)=0.
              PLM4(2)=0.
            ELSE IF (LTYP.EQ.'PO') THEN
              PRINT * , ' '
              PRINT * , 'Current latitude of 1st edge point: ',PLM1(1)
              PRINT * , 'Current longitude of 1st edge point:',PLM1(2)
              PRINT * , 'Current latitude of 2nd edge point: ',PLM2(1)
              PRINT * , 'Current longitude of 2nd edge point:',PLM2(2)
              PRINT * , 'Current latitude of 3rd edge point: ',PLM3(1)
              PRINT * , 'Current longitude of 3rd edge point:',PLM3(2)
              PRINT * , 'Current latitude of 4th edge point: ',PLM4(1)
              PRINT * , 'Current longitude of 4th edge point:',PLM4(2)
              PRINT * , ' '
              PRINT * , 'Enter new latitude of 1st edge point:'
              CALL EMRDRN (PLM1(1),PLM1(1))
              PRINT * , ' '
              PRINT * , 'Enter new longitude of 1st edge point:'
              CALL EMRDRN (PLM1(2),PLM1(2))
              PRINT * , ' '
              PRINT * , 'Enter new latitude of 2nd edge point:'
              CALL EMRDRN (PLM2(1),PLM2(1))
              PRINT * , ' '
              PRINT * , 'Enter new longitude of 2nd edge point:'
              CALL EMRDRN (PLM2(2),PLM2(2))
              PRINT * , ' '
              PRINT * , 'Enter new latitude of 3rd edge point:'
              CALL EMRDRN (PLM3(1),PLM3(1))
              PRINT * , ' '
              PRINT * , 'Enter new longitude of 3rd edge point:'
              CALL EMRDRN (PLM3(2),PLM3(2))
              PRINT * , ' '
              PRINT * , 'Enter new latitude of 4th edge point:'
              CALL EMRDRN (PLM4(1),PLM4(1))
              PRINT * , ' '
              PRINT * , 'Enter new longitude of 4th edge point:'
              CALL EMRDRN (PLM4(2),PLM4(2))
            ELSE IF (LTYP.EQ.'AN') THEN
              PRINT * , ' '
              PRINT * , 'Current angular distance to left:  ',PLM1(1)
              PRINT * , 'Current angular distance to right: ',PLM2(1)
              PRINT * , 'Current angular distance to bottom:',PLM3(1)
              PRINT * , 'Current angular distance to top:   ',PLM4(1)
              PRINT * , ' '
              PRINT * , 'Enter new angular distance to left:'
              CALL EMRDRN (PLM1(1),PLM1(1))
              PRINT * , ' '
              PRINT * , 'Enter new angular distance to right:'
              CALL EMRDRN (PLM2(1),PLM2(1))
              PRINT * , ' '
              PRINT * , 'Enter new angular distance to bottom:'
              CALL EMRDRN (PLM3(1),PLM3(1))
              PRINT * , ' '
              PRINT * , 'Enter new angular distance to top:'
              CALL EMRDRN (PLM4(1),PLM4(1))
              PLM1(2)=0.
              PLM2(2)=0.
              PLM3(2)=0.
              PLM4(2)=0.
            ELSE IF (LTYP.EQ.'LI') THEN
              PRINT * , ' '
              PRINT * , 'Current U coordinate of left edge:  ',PLM1(1)
              PRINT * , 'Current U coordinate of right edge: ',PLM2(1)
              PRINT * , 'Current V coordinate of bottom edge:',PLM3(1)
              PRINT * , 'Current V coordinate of top edge:   ',PLM4(1)
              PRINT * , ' '
              IF (PTYP.EQ.'LC'.OR.PTYP.EQ.'OR'.OR.PTYP.EQ.'SV') THEN
                PRINT * , 'U and V values should range from -1 to +1.'
              ELSE IF (PTYP.EQ.'ST'.OR.PTYP.EQ.'LE'.OR.PTYP.EQ.'GN')THEN
                PRINT * , 'U and V values should range from -2 to +2.'
              ELSE IF (PTYP.EQ.'AE'.OR.PTYP.EQ.'ME') THEN
                PRINT * , 'U and V values should range from -pi to +pi.'
              ELSE IF (PTYP.EQ.'CE') THEN
                PRINT * , 'U values should range from -180 to +180, V va
     +lues from -90 to +90.'
              ELSE IF (PTYP.EQ.'EA') THEN
                PRINT * , 'U values should range from -pi to +pi, V valu
     +es from -4/3 to +4/3.'
              ELSE IF (PTYP.EQ.'MO') THEN
                PRINT * , 'U values should range from -2 to +2, V values
     + from -1 to +1.'
              ELSE IF (PTYP.EQ.'RO') THEN
                PRINT * , 'U values should range from -1 to +1, V values
     + from -.5072 to +.5072.'
              ELSE IF (PTYP.EQ.'UT') THEN
                PRINT * , ' '
                PRINT * , 'The ranges of U and V values will depend on'
                PRINT * , 'the particular USGS transformation chosen;'
                PRINT * , 'UMIN, UMAX, VMIN, and VMAX from the last'
                PRINT * , 'EZMAP initialization were as follows:'
                PRINT * , ' '
                PRINT * , ' ',XWDL,XWDR,YWDB,YWDT
              END IF
              PRINT * , ' '
              PRINT * , 'Enter new U coordinate of left edge:'
              CALL EMRDRN (PLM1(1),PLM1(1))
              PRINT * , ' '
              PRINT * , 'Enter new U coordinate of right edge:'
              CALL EMRDRN (PLM2(1),PLM2(1))
              PRINT * , ' '
              PRINT * , 'Enter new V coordinate of bottom edge:'
              CALL EMRDRN (PLM3(1),PLM3(1))
              PRINT * , ' '
              PRINT * , 'Enter new V coordinate of top edge:'
              CALL EMRDRN (PLM4(1),PLM4(1))
              PLM1(2)=0.
              PLM2(2)=0.
              PLM3(2)=0.
              PLM4(2)=0.
            ELSE IF (LTYP.EQ.'GR') THEN
              PRINT * , ' '
              PRINT * , 'Latitude now at bottom of range: ',PLM1(1)
              PRINT * , 'Longitude now at bottom of range:',PLM2(1)
              PRINT * , 'Latitude now at top of range:    ',PLM3(1)
              PRINT * , 'Longitude now at top of range:   ',PLM4(1)
              PRINT * , ' '
              PRINT * , 'Enter new latitude at bottom of range:'
              CALL EMRDRN (PLM1(1),PLM1(1))
              PRINT * , ' '
              PRINT * , 'Enter new longitude at bottom of range:'
              CALL EMRDRN (PLM2(1),PLM2(1))
              PRINT * , ' '
              PRINT * , 'Enter new latitude at top of range:'
              CALL EMRDRN (PLM3(1),PLM3(1))
              PRINT * , ' '
              PRINT * , 'Enter new longitude at top of range:'
              CALL EMRDRN (PLM4(1),PLM4(1))
              PLM1(2)=0.
              PLM2(2)=0.
              PLM3(2)=0.
              PLM4(2)=0.
            END IF
          END IF
C
          PRINT * , ' '
          PRINT * , 'Change outline dataset parameters (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
            PRINT * , 'Possible outline dataset selectors:'
            PRINT * , ' '
            PRINT * , '  NO => None'
            PRINT * , '  CO => Continental outlines only'
            PRINT * , '  US => US outlines only'
            PRINT * , '  PS => Continental, political, and US outlines'
            PRINT * , '  PO => Continental and political outlines'
            PRINT * , '  E1 => The new database "Earth..1"'
            PRINT * , '  E2 => The new database "Earth..2"'
            PRINT * , '  E3 => The new database "Earth..3"'
            PRINT * , '  RG => The RANGS/GSHHS database'
            PRINT * , ' '
            PRINT * , 'Current outline dataset selector: ',ODNM
            PRINT * , ' '
            PRINT * , 'Enter new outline dataset selector:'
            READ  '(A2)', CTMP
            CALL MUPPER (CTMP)
            IF (CTMP.EQ.'NO'.OR.CTMP.EQ.'CO'.OR.CTMP.EQ.'US'.OR.
     +          CTMP.EQ.'PS'.OR.CTMP.EQ.'PO'.OR.CTMP.EQ.'E1'.OR.
     +          CTMP.EQ.'E2'.OR.CTMP.EQ.'E3'.OR.CTMP.EQ.'RG') ODNM=CTMP
            IF (ODNM.EQ.'E1'.OR.ODNM.EQ.'E2'.OR.ODNM.EQ.'E3') THEN
              PRINT * , ' '
              PRINT * , 'Current database level specifier:',ILVL
              PRINT * , ' '
              PRINT * , 'Enter new database level specifier:'
              CALL EMRDIN (ILVL,ILVL)
            ELSE IF (ODNM.EQ.'RG') THEN
              PRINT * , ' '
              PRINT * , 'Current name of data directory: ',USNM
              PRINT * , ' '
              PRINT * , 'Enter new name of data directory:'
              READ (*,'(A60)') CHRT
              IF (CHRT.NE.' ') USNM=CHRT
              PRINT * , ' '
              PRINT * , 'Current database level specifier:',MIN(5,IRGL)
              PRINT * , ' '
              PRINT * , 'Enter new database level specifier (0 => use'
              PRINT * , 'highest resolution, 1 to 4 => use successively'
              PRINT * , 'lower resolutions, 5 => pick a value for me):'
              CALL EMRDIN (IRGL,IRGL)
              IRGL=MAX(0,MIN(5,IRGL))
              PRINT * , ' '
              PRINT * , 'Change any color-fill flags (Y or N)?'
              READ  '(A1)', COMD
              IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
                PRINT * , ' '
                IF (LCSF(1).GE.0) THEN
                  PRINT * , 'Ocean will be filled.'
                ELSE
                  PRINT * , 'Ocean will not be filled.'
                END IF
                PRINT * , ' '
                PRINT * , 'Toggle this fill flag (Y or N)?'
                READ  '(A1)', COMD
                IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') LCSF(1)=-LCSF(1)
                PRINT * , ' '
                IF (LCSF(2).GE.0) THEN
                  PRINT * , 'Land will be filled.'
                ELSE
                  PRINT * , 'Land will not be filled.'
                END IF
                PRINT * , ' '
                PRINT * , 'Toggle this fill flag (Y or N)?'
                READ  '(A1)', COMD
                IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') LCSF(2)=-LCSF(2)
                PRINT * , ' '
                IF (LCSF(3).GE.0) THEN
                  PRINT * , 'Lakes will be filled.'
                ELSE
                  PRINT * , 'Lakes will not be filled.'
                END IF
                PRINT * , ' '
                PRINT * , 'Toggle this fill flag (Y or N)?'
                READ  '(A1)', COMD
                IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') LCSF(3)=-LCSF(3)
                PRINT * , ' '
                IF (LCSF(4).GE.0) THEN
                  PRINT * , 'Islands in lakes will be filled.'
                ELSE
                  PRINT * , 'Islands in lakes will not be filled.'
                END IF
                PRINT * , ' '
                PRINT * , 'Toggle this fill flag (Y or N)?'
                READ  '(A1)', COMD
                IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') LCSF(4)=-LCSF(4)
                PRINT * , ' '
                IF (LCSF(5).GE.0) THEN
                  PRINT * , 'Ponds on islands in lakes will be filled.'
                ELSE
                  PRINT * , 'Ponds on islands in lakes will not be fille
     +d.'
                END IF
                PRINT * , ' '
                PRINT * , 'Toggle this fill flag (Y or N)?'
                READ  '(A1)', COMD
                IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') LCSF(5)=-LCSF(5)
                CALL MDRGSC (LCOL,LCSF)
              END IF
            END IF
          END IF
C
C Other parameters to provide access to:
C
C   GD LS
C
          PRINT * , ' '
          PRINT * , 'Change miscellaneous minor EZMAP parameters (Y or N
     +)?'
          READ  '(A1)', COMD
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
            PRINT * , ' '
            PRINT * , 'Current perimeter flag:           ',IPER
            PRINT * , 'Current elliptical-perimeter flag:',IELL
            PRINT * , 'Current labelling flag:           ',ILBL
            PRINT * , 'Current dotted-outline flag:      ',IDOT
            PRINT * , 'Current distance between dots:    ',IDBD
            PRINT * , 'Current grid spacing:             ',GSPC
            PRINT * , 'Current latitude line spacing:    ',GSLA
            PRINT * , 'Current longitude line spacing:   ',GSLO
            PRINT * , 'Current polar latitude control:   ',GLAT
            PRINT * , 'Current polar longitude control:  ',GLON
            PRINT * , 'Current RG interpolation control: ',IINT
            PRINT * , 'Current RG data processing flag:  ',IRPF
            PRINT * , ' '
            PRINT * , 'Enter new perimeter flag (0 or 1):'
            CALL EMRDIN (IPER,IPER)
            IPER=MAX(0,MIN(1,IPER))
            PRINT * , ' '
            PRINT * , 'Enter new elliptical-perimeter flag (0 or 1):'
            CALL EMRDIN (IELL,IELL)
            IELL=MAX(0,MIN(1,IELL))
            PRINT * , ' '
            PRINT * , 'Enter new labelling flag (0 or 1):'
            CALL EMRDIN (ILBL,ILBL)
            ILBL=MAX(0,MIN(1,ILBL))
            PRINT * , ' '
            PRINT * , 'Enter new dotted-outline flag (0 or 1):'
            CALL EMRDIN (IDOT,IDOT)
            IDOT=MAX(0,MIN(1,IDOT))
            PRINT * , ' '
            PRINT * , 'Enter new distance between dots (out of 32768):'
            CALL EMRDIN (IDBD,IDBD)
            IDBD=MAX(1,IDBD)
            PRINT * , ' '
            PRINT * , 'Enter new grid spacing (0 turns grid off):'
            CALL EMRDRN (GSPC,GSPC)
            IGRD=MAX(0.,MIN(90.,GSPC))
            PRINT * , ' '
            PRINT * ,  'Enter new latitude line spacing (zero if same as
     + grid spacing):'
            CALL EMRDRN (GSLA,GSLA)
            GSLA=MAX(0.,MIN(90.,GSLA))
            PRINT * , ' '
            PRINT * , 'Enter new longitude line spacing (zero if same as
     + grid spacing):'
            CALL EMRDRN (GSLO,GSLO)
            GSLO=MAX(0.,MIN(90.,GSLO))
            PRINT * , ' '
            PRINT * , 'Enter latitude at edge of polar region (integer f
     +rom 0 to 90):'
            CALL EMRDRN (GLAT,GLAT)
            GLAT=MAX(0.,MIN(90.,REAL(INT(GLAT))))
            PRINT * , ' '
            PRINT * , 'Enter polar-region longitude-line spacing (0 to s
     +uppress them at poles):'
            CALL EMRDRN (GLON,GLON)
            GLON=MAX(0.,MIN(90.,GLON))
            PRINT * , ' '
            PRINT * , 'Enter RANGS/GSHHS interpolation flag:'
            CALL EMRDIN (IINT,IINT)
            IINT=MAX(1,MIN(999,IINT))
            PRINT * , ' '
            PRINT * , 'Enter RANGS/GSHHS data processing flag:'
            CALL EMRDIN (IRPF,IRPF)
            IRPF=MAX(0,MIN(2,IRPF))
          END IF
C
          PRINT * , ' '
          PRINT * , 'Change miscellaneous (non-EZMAP) parameters'
          PRINT * , 'within this demo program (Y or N)?'
          READ  '(A1)', COMD
C
          IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
C
            PRINT * , ' '
            IF (IBOW.EQ.0) THEN
              PRINT * , 'The color scheme flag is 0 (white on black).'
            ELSE
              PRINT * , 'The color scheme flag is 1 (black on white).'
            END IF
            PRINT * , ' '
            PRINT * , 'Enter new color scheme flag (0=>WOB, 1=>BOW):'
            CALL EMRDIN (ITMP,IBOW)
            ITMP=MAX(0,MIN(1,ITMP))
            IF (ITMP.NE.IBOW) THEN
              IBOW=ITMP
              CALL GCLRWK (1,1)
              DO 107 IWID=1,3
                IF (IWID.NE.1) CALL GACWK (IWID)
                IF (IBOW.EQ.0) THEN
                  CALL GSCR (IWID,  0,0.,0.,0.)
                  CALL GSCR (IWID,  1,1.,1.,1.)
                ELSE
                  CALL GSCR (IWID,  0,1.,1.,1.)
                  CALL GSCR (IWID,  1,0.,0.,0.)
                END IF
                CALL GSCR (IWID,  2,1.,0.,0.)
                CALL GSCR (IWID,  3,0.,1.,0.)
                CALL GSCR (IWID,  4,0.,0.,1.)
                CALL GSCR (IWID,  5,0.,1.,1.)
                CALL GSCR (IWID,  6,1.,0.,1.)
                CALL GSCR (IWID,  7,1.,1.,0.)
                CALL GSCR (IWID,  8,.7,.7,0.)
                CALL GSCR (IWID,  9,1.,1.,0.)
                CALL GSCR (IWID, 10,0.,0.,0.)
                CALL GSCR (IWID, 11,1.,1.,.6)
                CALL GSCR (IWID, 12,1.,0.,0.)
                CALL GSCR (IWID,101,.2,.2,.8)
                CALL GSCR (IWID,102,.6,.4,.2)
                CALL GSCR (IWID,103,.2,.6,.4)
                CALL GSCR (IWID,104,.2,.8,.2)
                CALL GSCR (IWID,105,.4,.6,.2)
                CALL GSCR (IWID,106,.2,.4,.6)
                CALL GSCR (IWID,107,.6,.6,.6)
                IF (IWID.NE.1) CALL GDAWK (IWID)
  107         CONTINUE
            END IF
C
            PRINT * , ' '
            IF (ISTY.EQ.0) THEN
              PRINT * , 'The style flag is 0 (simple map).'
            ELSE
              PRINT * , 'The style flag is 1 (color-filled map).'
            END IF
            PRINT * , ' '
            PRINT * , 'Enter new style flag (0 => simple map, 1 => color
     +-filled map):'
            CALL EMRDIN (ISTY,ISTY)
            ISTY=MAX(0,MIN(1,ISTY))
C
            PRINT * , ' '
            IF (IUXW.EQ.0) THEN
              PRINT * , 'The X-workstation-update flag is 0 (off).'
            ELSE
              PRINT * , 'The X-workstation-update flag is 1 (on).'
            END IF
            PRINT * , ' '
            PRINT * , 'Enter new X-workstation-update flag (0=>off, 1=>o
     +n):'
            CALL EMRDIN (IUXW,IUXW)
            IUXW=MAX(0,MIN(1,IUXW))
C
            PRINT * , ' '
            IF (ILPT.EQ.0) THEN
              PRINT * , 'The informational label flag is 0 (no label).
     +'
            ELSE
              PRINT * , 'The informational label flag is 1 (do label).
     +'
            END IF
            PRINT * , ' '
            PRINT * , 'Enter new informational label flag (0=>off, 1=>on
     +):'
            CALL EMRDIN (ILPT,ILPT)
            ILPT=MAX(0,MIN(1,ILPT))
C
            PRINT * , ' '
            IF (LABL.EQ.' ') THEN
              PRINT * , 'No plot label is currently defined.'
            ELSE
              PRINT * , 'The following plot label is currently defined:'
              WRITE (*,'(2X,A60)') LABL
            END IF
            PRINT * , ' '
            PRINT * , 'Enter one of the following:'
            PRINT * , ' '
            PRINT * , '  A blank line (use the current label)'
            PRINT * , '  A new label (60 or fewer characters)'
            PRINT * , '  The word "NONE" or "none" (no label)'
            PRINT * , ' '
            READ (*,'(A60)') CHRT
            IF (CHRT.EQ.'NONE'.OR.CHRT.EQ.'none') THEN
              LABL=' '
            ELSE IF (CHRT.NE.' ') THEN
              LABL=CHRT
            END IF
C
            PRINT * , ' '
            IF (ILLB.EQ.0) THEN
              PRINT * , 'The lat/lon label flag is 0 (no labels).'
            ELSE IF (ILLB.EQ.-1) THEN
              PRINT * , 'The lat/lon label flag is -1 (labels outside ma
     +p edge).'
            ELSE IF (ILLB.EQ.-1) THEN
              PRINT * , 'The lat/lon label flag is +1 (labels inside map
     + edge.)'
            END IF
            PRINT * , 'Current width of characters in labels: ',CSLL
            PRINT * , ' '
            PRINT * , 'Enter new lat/lon label flag:'
            PRINT * , ' '
            PRINT * , '   0=>no lat/lon labels'
            PRINT * , '  -1=>labels outside map edge'
            PRINT * , '  +1=>labels inside map edge'
            PRINT * , ' '
            CALL EMRDIN (ILLB,ILLB)
            ILLB=MAX(-1,MIN(+1,ILLB))
            PRINT * , ' '
            PRINT * , 'Enter new character width:'
            PRINT * , ' '
            CALL EMRDRN (CSLL,CSLL)
            CSLL=MAX(.001,MIN(.1,CSLL))
C
            PRINT * , ' '
            IF (ISPG.LT.0) THEN
              PRINT * , 'State Plane zones will not be drawn.'
            ELSE IF (ISPG.EQ.0) THEN
              PRINT * , 'State Plane zones for NAD1927 will be drawn.'
            ELSE
              PRINT * , 'State Plane zones for NAD1983 will be drawn.'
            END IF
            PRINT * , ' '
            PRINT * , 'Enter new State Plane zone flag:'
            PRINT * , ' '
            PRINT * , '  -n => off'
            PRINT * , '   0 => NAD1927'
            PRINT * , '  +n => NAD1983'
            CALL EMRDIN (ISPG,ISPG)
            IF (ISPG.LT.0) THEN
              ISPG=-1
            ELSE IF (ISPG.GT.0) THEN
              ISPG=8
            END IF
            IF (ISPG.GE.0) THEN
              PRINT * , ' '
              PRINT * , 'Current value of character width for labels:',
     +CHSZ
              PRINT * , ' '
              PRINT * , 'Enter new character width (from .01 to 10 degre
     +es):'
              CALL EMRDRN (CHSZ,CHSZ)
              CHSZ=MAX(.01,MIN(10.,CHSZ))
            END IF
C
          END IF
C
          GO TO 105
C
        ELSE IF (COMD.EQ.'P'.OR.COMD.EQ.'p') THEN
C
          PRINT * , ' '
          PRINT * , 'Selecting a new predefined example plot.  Enter'
          PRINT * , 'one of the following:'
          PRINT * , ' '
          PRINT * , '  N => next plot in the list'
          PRINT * , '  P => previous plot in the list'
          PRINT * , '  S => specified plot in the list'
          PRINT * , ' '
          PRINT * , '  E => expand current plot into a square'
C
          READ  '(A1)', COMD
C
          IF      (COMD.EQ.'E'.OR.COMD.EQ.'e') THEN
            PRINT * , ' '
            PRINT * , 'Expanding the current plot into a square.'
            LTYP='LI'
            PLM1(1)=.5*(XWDL+XWDR)-.5*MAX(XWDR-XWDL,YWDT-YWDB)
            PLM2(1)=.5*(XWDL+XWDR)+.5*MAX(XWDR-XWDL,YWDT-YWDB)
            PLM3(1)=.5*(YWDB+YWDT)-.5*MAX(XWDR-XWDL,YWDT-YWDB)
            PLM4(1)=.5*(YWDB+YWDT)+.5*MAX(XWDR-XWDL,YWDT-YWDB)
            PLM1(2)=0.
            PLM2(2)=0.
            PLM3(2)=0.
            PLM4(2)=0.
            GO TO 105
          ELSE IF (COMD.EQ.'N'.OR.COMD.EQ.'n') THEN
            PRINT * , ' '
            PRINT * , 'Selected next one in list.'
          ELSE IF (COMD.EQ.'P'.OR.COMD.EQ.'p') THEN
            PRINT * , ' '
            PRINT * , 'Selected previous one in list.'
            IPDP=MOD(IPDP+65,68)+1
          ELSE IF (COMD.EQ.'S'.OR.COMD.EQ.'s') THEN
            PRINT * , ' '
            PRINT * , 'Enter an integer between 1 and 68 (use 1-13 for'
            PRINT * , 'original EZMAP projection types, 14-46 for USGS'
            PRINT * , 'transformations, and 47-68 for State Plane Grid'
            PRINT * , 'plots):'
            PRINT * , ' '
            PRINT * , '   1 => LC     24 => UT-06  47 => SPG-US-NE--27'
            PRINT * , '   2 => GN     25 => UT-06  48 => SPG-US-NE--83'
            PRINT * , '   3 => ME     26 => UT-07  49 => SPG-US-SE--27'
            PRINT * , '   4 => ST     27 => UT-08  50 => SPG-US-SE--83'
            PRINT * , '   5 => LE     28 => UT-08  51 => SPG-US-NEC-27'
            PRINT * , '   6 => OR     29 => UT-08  52 => SPG-US-NEC-83'
            PRINT * , '   7 => SV     30 => UT-08  53 => SPG-US-SEC-27'
            PRINT * , '   8 => AE     31 => UT-09  54 => SPG-US-SEC-83'
            PRINT * , '   9 => RO     32 => UT-10  55 => SPG-US-NWC-27'
            PRINT * , '  10 => MO     33 => UT-11  56 => SPG-US-NWC-83'
            PRINT * , '  11 => CE     34 => UT-12  57 => SPG-US-SWC-27'
            PRINT * , '  12 => EA     35 => UT-13  58 => SPG-US-SWC-83'
            PRINT * , '  13 => RM     36 => UT-14  59 => SPG-US-NW--27'
            PRINT * , '  14 => UT-01  37 => UT-15  60 => SPG-US-NW--83'
            PRINT * , '  15 => UT-01  38 => UT-16  61 => SPG-US-SW--27'
            PRINT * , '  16 => UT-02  39 => UT-17  62 => SPG-US-SW--83'
            PRINT * , '  17 => UT-02  40 => UT-18  63 => SPG-ALASKA-27'
            PRINT * , '  18 => UT-02  41 => UT-19  64 => SPG-ALASKA-83'
            PRINT * , '  19 => UT-03  42 => UT-20  65 => SPG-HAWAII-27'
            PRINT * , '  20 => UT-03  43 => UT-20  66 => SPG-HAWAII-83'
            PRINT * , '  21 => UT-04  44 => UT-21  67 => SPG-ALL----27'
            PRINT * , '  22 => UT-04  45 => UT-22  68 => SPG-ALL----83'
            PRINT * , '  23 => UT-05  46 => UT-23                     '
            PRINT * , ' '
            CALL EMRDIN (IPDP,IPDP)
            IPDP=MAX(1,MIN(68,IPDP))
          END IF
C
          PRINT * , ' '
          PRINT * , 'Number of predefined example plot selected:',IPDP
C
          CALL INNPDP (IPDP,PTYP,PLAT,PLON,ROTA,LTYP,PLM1,PLM2,PLM3,
     +                 PLM4,IPRJ,IZON,ISPH,PARA,IDSL,ISPG,ODNM,ILVL,
     +                 CHSZ,LABL,ILLB)
C
          GO TO 105
C
        ELSE IF (COMD.EQ.'Q'.OR.COMD.EQ.'q') THEN
C
          CALL GDAWK (1)
          CALL GCLWK (1)
          CALL GCLWK (2)
          IF (NFPS.NE.0) THEN
            CALL GACWK (3)
            CALL GCLRWK (3,1)
            CALL GDAWK (3)
          END IF
          CALL GCLWK (3)
          CALL GCLKS
          STOP
C
        ELSE IF (COMD.EQ.'R'.OR.COMD.EQ.'r') THEN
C
          PRINT * , ' '
          PRINT * , 'Redrawing the picture.'
C
          GO TO 105
C
        ELSE IF (COMD.EQ.'S'.OR.COMD.EQ.'s') THEN
C
          PRINT * , ' '
          PRINT * , 'Save the current frame.  Note that parts of the'
          PRINT * , 'image resulting from "T" commands are not saved.'
          PRINT * , ' '
          PRINT * , 'One may save to an NCGM, to a Postscript file, or'
          PRINT * , 'both.'
          PRINT * , ' '
          PRINT * , 'Do you want to save to an NCAR CGM file (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.NE.'Y'.AND.COMD.NE.'y') THEN
            ISNG=0
          ELSE
            ISNG=1
            CALL GACWK (2)
            NFNG=NFNG+1
            IF (NFNG.GT.1) CALL GCLRWK (2,1)
          END IF
          PRINT * , ' '
          PRINT * , 'Do you want to save to a PostScript file (Y or N)?'
          READ  '(A1)', COMD
          IF (COMD.NE.'Y'.AND.COMD.NE.'y') THEN
            ISPS=0
          ELSE
            ISPS=1
            CALL GACWK (3)
            NFPS=NFPS+1
            IF (NFPS.GT.1) CALL GCLRWK (3,1)
          END IF
          IF (ISNG.NE.0.OR.ISPS.NE.0) THEN
            CALL GDAWK(1)
            CALL MAPPOS (XVPL,XVPR,YVPB,YVPT)
            CALL MAPROJ (PTYP,PLAT,PLON,ROTA)
            IF (PTYP.EQ.'SV') THEN
              CALL MPSETR ('SA',SALT)
              CALL MPSETR ('S1',SAN1)
              CALL MPSETR ('S2',SAN2)
            END IF
            CALL MAPSET (LTYP,PLM1,PLM2,PLM3,PLM4)
            CALL MPSETI ('PE',IPER)
            CALL MPSETI ('EL',IELL)
            CALL MPSETI ('LA',ILBL)
            CALL MPSETI ('DO',IDOT)
            CALL MPSETI ('DD',IDBD)
            CALL MPSETR ('GR',GSPC)
            CALL MPSETR ('GT',GSLA)
            CALL MPSETR ('GN',GSLO)
            CALL MPSETR ('GP',1000.*GLAT+GLON)
            CALL MPSETI ('II',1000*IINT+IINT)
            CALL MPSETI ('RP',IRPF)
            IF (PTYP.NE.'UT') THEN
              CALL MAPINT
            ELSE
              CALL MPUTIN (IPRJ,IZON,ISPH,PARA,0.D0,0.D0,0.D0,0.D0)
            END IF
            IF (ISTY.EQ.0) THEN
              CALL MAPGRD
              CALL MAPLBL
              IF (ODNM.EQ.'E1'.OR.ODNM.EQ.'E2'.OR.ODNM.EQ.'E3') THEN
                IF (ODNM.EQ.'E1') CALL MPLNDR ('Earth..1',ILVL)
                IF (ODNM.EQ.'E2') CALL MPLNDR ('Earth..2',ILVL)
                IF (ODNM.EQ.'E3') CALL MPLNDR ('Earth..3',ILVL)
              ELSE IF (ODNM.EQ.'RG') THEN
                CALL MDRGOL (MOD(IRGL,5 ),RWRK,LRWK)
              ELSE
                CALL MAPLOT
              END IF
              CALL DRSPGD (ISPG,IZN0,PLT0,PLN0,CLT0,CLN0,
     +                          IZN8,PLT8,PLN8,CLT8,CLN8,CHSZ)
            ELSE
              IF (ODNM.EQ.'E1'.OR.ODNM.EQ.'E2'.OR.ODNM.EQ.'E3') THEN
                CALL MPSETI ('VS',20)
                CALL ARINAM (IAMA,LAMA)
                IF (ODNM.EQ.'E1') CALL MPLNAM ('Earth..1',ILVL,IAMA)
                IF (ODNM.EQ.'E2') CALL MPLNAM ('Earth..2',ILVL,IAMA)
                IF (ODNM.EQ.'E3') CALL MPLNAM ('Earth..3',ILVL,IAMA)
                CALL ARSCAM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,COLORB)
                IF (ODNM.EQ.'E1') CALL MPLNDR ('Earth..1',ILVL)
                IF (ODNM.EQ.'E2') CALL MPLNDR ('Earth..2',ILVL)
                IF (ODNM.EQ.'E3') CALL MPLNDR ('Earth..3',ILVL)
                CALL MAPGRM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,DRAWLB)
              ELSE IF (ODNM.EQ.'RG') THEN
                CALL MDRGSF (MOD(IRGL,5),RWRK,LRWK,IAMA,LAMA)
                CALL MAPGRD
              ELSE
                CALL MPSETC ('OU',ODNM)
                CALL MPSETI ('VS',0)
                CALL ARINAM (IAMA,LAMA)
                CALL MAPBLA (IAMA)
                CALL ARSCAM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,COLORA)
                CALL MAPLOT
                CALL MAPGRM (IAMA,XCRA,YCRA,NCRA,IAAI,IAGI,NGPS,DRAWLA)
              END IF
              CALL MAPLMB
              CALL MAPLBL
              CALL DRSPGD (ISPG,IZN0,PLT0,PLN0,CLT0,CLN0,
     +                          IZN8,PLT8,PLN8,CLT8,CLN8,CHSZ)
            END IF
            IF (ILLB.NE.0) THEN
              P=REAL(SIGN(1,ILLB))
              CALL PCSETI ('BF - BOX FLAG',0)
              CALL MDLBLT(RVPL,RVPB,RVPL,RVPT,+P*CSLL,    0.,CSLL,0.,-P)
              CALL MDLBLT(RVPR,RVPB,RVPR,RVPT,-P*CSLL,    0.,CSLL,0.,+P)
              CALL MDLBLN(RVPL,RVPB,RVPR,RVPB,0.,+1.5*P*CSLL,CSLL,0.,0.)
              CALL MDLBLN(RVPL,RVPT,RVPR,RVPT,0.,-1.5*P*CSLL,CSLL,0.,0.)
              CALL PCSETI ('BF - BOX FLAG',3)
            END IF
            IF (LABL.NE.' ') CALL PLCHHQ (CFUX(.5),CFUY(.965),
     +                                    LABL(1:MPILNB(LABL)),
     +                                              .018,0.,0.)
            IF (ILPT.NE.0) CALL DRPTLB (PTYP,IPRJ,IZON,ISPH,PARA,IPRF,
     +                                                      ODNM,ILVL)
            CALL SFLUSH
            IF (ISNG.NE.0) CALL GDAWK (2)
            IF (ISPS.NE.0) CALL GDAWK (3)
            CALL GACWK(1)
          END IF
C
          GO TO 106
C
        ELSE IF (COMD.EQ.'T'.OR.COMD.EQ.'t') THEN
C
          PRINT * , ' '
          PRINT * , 'Test some aspect of EZMAP.'
          PRINT * , ' '
          PRINT * , 'Enter one of the following commands:'
          PRINT * , ' '
          PRINT * , '  C  =>  Contour test'
          PRINT * , '  I  =>  Inverse routine test'
          PRINT * , '  S  =>  Spiral test'
C
          READ  '(A1)', COMD
C
          IF (COMD.EQ.'C'.OR.COMD.EQ.'c') THEN
            PRINT * , ' '
            PRINT * , 'This test uses EZMAP''s inverse transformation'
            PRINT * , 'routines to generate large arrays of latitudes'
            PRINT * , 'and longitudes as functions of the fractional'
            PRINT * , 'coordinates of points in the plotter frame.  The'
            PRINT * , 'arrays are contoured.  The contour lines should'
            PRINT * , 'be observed to lie on top of the lat/lon grid'
            PRINT * , 'previously drawn.'
            PRINT * , ' '
            PRINT * , 'As the latitudes and longitudes are computed,'
            PRINT * , 'points where the inverse transformation is not'
            PRINT * , 'defined are marked with little X''s whose color'
            PRINT * , 'implies the reason the inverse was not defined.'
            PRINT * , ' '
            PRINT * , 'One may change the color of the contours and the'
            PRINT * , 'dimensions of the arrays used.  Using larger'
            PRINT * , 'arrays will more accurately duplicate the grid'
            PRINT * , 'line positions, but is time-consuming.'
            PRINT * , ' '
            PRINT * , 'Change contouring parameters (Y or N)?:'
            READ  '(A1)', COMD
            IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
              PRINT * , ' '
              PRINT * , 'Current contouring color:          ',ICNC
              PRINT * , 'Current dimension of contour array:',NCDM
              PRINT * , ' '
              PRINT * , 'Enter new contouring color (between 0 and 7):'
              PRINT * , ' '
              IF (IBOW.EQ.0) THEN
                PRINT * , '  0 => black'
                PRINT * , '  1 => white'
              ELSE
                PRINT * , '  0 => white'
                PRINT * , '  1 => black'
              END IF
              PRINT * , '  2 => red'
              PRINT * , '  3 => green'
              PRINT * , '  4 => blue'
              PRINT * , '  5 => cyan'
              PRINT * , '  6 => magenta'
              PRINT * , '  7 => yellow'
              CALL EMRDIN (ICNC,ICNC)
              ICNC=MAX(0,MIN(7,ICNC))
              PRINT * , ' '
              PRINT * , 'Enter new contour array dimension (between 51 a
     +nd 501):'
              CALL EMRDIN (NCDM,NCDM)
              NCDM=MAX(51,MIN(501,NCDM))
            END IF
            PRINT * , ' '
            PRINT * , 'Contouring latitudes and longitudes.'
            CALL CONINV (ICNC,NCDM)
            CALL MAPINT
            GO TO 106
          ELSE IF (COMD.EQ.'I'.OR.COMD.EQ.'i') THEN
            IF (ISTY.NE.0.AND.ODNM.NE.'RG') THEN
              CALL GSPLCI (2)
              DO 109 I=0,9
                CALL PLCHMQ (CFUX(REAL(I)/10.+.05),CFUY(.05),
     +                       CHAR(ICHAR('0')+I),.015,0.,0.)
                CALL PLCHMQ (CFUX(.05),CFUY(REAL(I)/10.+.05),
     +                       CHAR(ICHAR('0')+I),.015,0.,0.)
  109         CONTINUE
              PRINT * , ' '
              PRINT * , 'This test allows one to repeatedly select one'
              PRINT * , 'of 100 subsquares of the plotter frame; each'
              PRINT * , 'subsquare so selected is further subdivided'
              PRINT * , 'into 100 subsubsquares and each of those is'
              PRINT * , 'filled with a color.  The color is selected'
              PRINT * , 'by a process that thoroughly tests the EZMAP'
              PRINT * , 'transformation routines (in both directions).'
              PRINT * , 'The effect should be to leave the map colored'
              PRINT * , 'as it was by the initial color fill.'
  110         PRINT * , ' '
              PRINT * , 'Enter X position (0-9, -1 to quit):'
              CALL EMRDIN (IHOR,0)
              IHOR=MAX(-1,MIN(9,IHOR))
              IF (IHOR.EQ.-1) GO TO 106
              PRINT * , ' '
              PRINT * , 'Enter Y position (0-9, -1 to quit):'
              CALL EMRDIN (IVER,0)
              IVER=MAX(-1,MIN(9,IVER))
              IF (IVER.EQ.-1) GO TO 106
              ICSF=1
              DO 113 I=0,99
                DO 112 J=0,99
                  XBOX(1)=CFUX(REAL(IHOR)/10.+REAL(I  )/1000.)
                  YBOX(1)=CFUY(REAL(IVER)/10.+REAL(J  )/1000.)
                  XBOX(2)=CFUX(REAL(IHOR)/10.+REAL(I+1)/1000.)
                  YBOX(2)=CFUY(REAL(IVER)/10.+REAL(J  )/1000.)
                  XBOX(3)=CFUX(REAL(IHOR)/10.+REAL(I+1)/1000.)
                  YBOX(3)=CFUY(REAL(IVER)/10.+REAL(J+1)/1000.)
                  XBOX(4)=CFUX(REAL(IHOR)/10.+REAL(I  )/1000.)
                  YBOX(4)=CFUY(REAL(IVER)/10.+REAL(J+1)/1000.)
                  XPOS=CFUX(REAL(IHOR)/10.+REAL(I  )/1000.+.0005)
                  YPOS=CFUY(REAL(IVER)/10.+REAL(J+1)/1000.+.0005)
                  CALL MAPTRI (XPOS,YPOS,RLAT,RLON)
                  IF (RLAT.NE.1.E12.AND.RLON.NE.1.E12) THEN
                    CALL MAPTRA (RLAT,RLON,XPOS,YPOS)
                    IF (XPOS.NE.1.E12.AND.YPOS.NE.1.E12) THEN
                      CALL ARGTAI (IAMA,XPOS,YPOS,IAAI,IAGI,NGPS,NAIR,
     +                                                           ICSF)
                      ICSF=0
                      IAI1=-1
                      DO 111 IAIR=1,NAIR
                        IF (IAGI(IAIR).EQ.1) IAI1=IAAI(IAIR)
  111                 CONTINUE
                      IF (IAI1.GT.0) THEN
                        IF (ODNM.NE.'E1'.AND.ODNM.NE.'E2'.AND.
     +                                       ODNM.NE.'E3') THEN
                          CALL GSFACI (MAPACI(IAI1)+100)
                        ELSE
                          CALL GSFACI (MPISCI(IAI1)+100)
                        END IF
                        CALL GFA (4,XBOX,YBOX)
                      END IF
                    END IF
                  END IF
  112           CONTINUE
  113         CONTINUE
              CALL PLOTIF (REAL(IHOR  )/10.,REAL(IVER  )/10.,0)
              CALL PLOTIF (REAL(IHOR+1)/10.,REAL(IVER  )/10.,1)
              CALL PLOTIF (REAL(IHOR+1)/10.,REAL(IVER+1)/10.,1)
              CALL PLOTIF (REAL(IHOR  )/10.,REAL(IVER+1)/10.,1)
              CALL PLOTIF (REAL(IHOR  )/10.,REAL(IVER  )/10.,1)
              CALL SFLUSH
              CALL GUWK (1,0)
              GO TO 110
            ELSE
              PRINT * , ' '
              PRINT * , 'Can''t do that test with "style" flag = 0 or'
              PRINT * , 'when using outline dataset selector "RG".'
              GO TO 106
            END IF
          ELSE IF (COMD.EQ.'S'.OR.COMD.EQ.'s') THEN
            PRINT * , ' '
            PRINT * , 'This test allows one to draw the projection of'
            PRINT * , 'a spiral beginning at a specified point on the'
            PRINT * , 'globe and ending at the point exactly opposite'
            PRINT * , 'that initial point.  One may change the color of'
            PRINT * , 'the spiral, the number of complete turns that it'
            PRINT * , 'makes, and the location of its starting point.'
            PRINT * , 'This tests the forward-transformation routines'
            PRINT * , 'and the proper drawing of lines as they appear'
            PRINT * , 'and disappear at the edges of a map.'
            PRINT * , ' '
            PRINT * , 'Change parameters of the spiral (Y or N)?:'
            READ  '(A1)', COMD
            IF (COMD.EQ.'Y'.OR.COMD.EQ.'y') THEN
              PRINT * , ' '
              PRINT * , 'Current spiral color:      ',ISPC
              PRINT * , 'Current number of turns:   ',NSPT
              PRINT * , 'Current latitude at start: ',SPLT
              PRINT * , 'Current longitude at start:',SPLN
              PRINT * , ' '
              PRINT * , 'Enter new spiral color (between 0 and 7):'
              PRINT * , ' '
              IF (IBOW.EQ.0) THEN
                PRINT * , '  0 => black'
                PRINT * , '  1 => white'
              ELSE
                PRINT * , '  0 => white'
                PRINT * , '  1 => black'
              END IF
              PRINT * , '  2 => red'
              PRINT * , '  3 => green'
              PRINT * , '  4 => blue'
              PRINT * , '  5 => cyan'
              PRINT * , '  6 => magenta'
              PRINT * , '  7 => yellow'
              CALL EMRDIN (ISPC,ISPC)
              ISPC=MAX(0,MIN(7,ISPC))
              PRINT * , ' '
              PRINT * , 'Enter new number of turns (between 1 and 1000):
     +'
              CALL EMRDIN (NSPT,NSPT)
              NSPT=MAX(1,MIN(1000,NSPT))
              PRINT * , ' '
              PRINT * , 'Enter new latitude at start (between -90 and 90
     +):'
              CALL EMRDRN (SPLT,SPLT)
              SPLT=MAX(-90.,MIN(90.,SPLT))
              PRINT * , ' '
              PRINT * ,  'Enter new longitude at start (between -180 and
     + 180):'
              CALL EMRDRN (SPLN,SPLN)
              SPLN=MAX(-180.,MIN(180.,SPLN))
            END IF
            PRINT * , ' '
            PRINT * , 'Drawing spiral.'
            CALL DRAWSP (ISPC,SPLT,SPLN,NSPT)
            GO TO 106
          ELSE
            PRINT * , ' '
            PRINT * , 'Illegal test command.'
            GO TO 106
          END IF
C
        ELSE
C
          PRINT * , ' '
          PRINT * , 'Illegal command.  Try again.'
          GO TO 106
C
        END IF
C
C An error occurred in a graphics call.  Inform user and go back to
C let him enter another command.
C
  901   PRINT * , ' '
        PRINT * , 'An error has occurred in a graphics call.  The error'
        PRINT * , 'message was as follows:'
        PRINT * , ' '
        CALL EPRIN
        PRINT * , ' '
        PRINT * , 'You will probably have to change some parameter to'
        PRINT * , 'avoid this error.'
        CALL ERROF
        GO TO 106
C
      END



      SUBROUTINE EMRDIN (IVAL,IDEF)
C
C This routine reads a line of input to get an integer value IVAL.  If
C the line is entirely blank, the default value IDEF is returned.  If
C the line contains an illegal character, another line is obtained.
C
        CHARACTER*128 LINE
C
C Assume the default value will be returned until we find otherwise.
C
        IVAL=IDEF
C
C Get a line from standard input.
C
        READ '(A128)' , LINE
C
        IF (LINE.NE.' ') THEN
          READ (LINE,*) IVAL
        END IF
C
C Done.
C
        RETURN
C
      END



      SUBROUTINE EMRDRN (RVAL,RDEF)
C
C This routine reads a line of input to get a real number RVAL.  If
C the line is entirely blank, the default value RDEF is returned.  If
C the line contains an illegal character, another line is obtained.
C
        CHARACTER*128 LINE
C
C Assume the default value will be returned until we find otherwise.
C
        RVAL=RDEF
C
C Get a line from standard input.
C
        READ '(A128)' , LINE
C
        IF (LINE.NE.' ') THEN
          READ (LINE,*) RVAL
        END IF
C
C Done.
C
        RETURN
C
      END



      SUBROUTINE COLORA (XCRA,YCRA,NCRA,IAAI,IAGI,NGPS)
        DIMENSION XCRA(NCRA),YCRA(NCRA),IAAI(NGPS),IAGI(NGPS)
        IAI1=-1
        DO 101 IGRP=1,NGPS
          IF (IAGI(IGRP).EQ.1) IAI1=IAAI(IGRP)
  101   CONTINUE
        IF (IAI1.GT.0) THEN
          CALL GSFACI (MAPACI(IAI1)+100)
          CALL GFA    (NCRA-1,XCRA,YCRA)
        END IF
        RETURN
      END



      SUBROUTINE DRAWLA (XCRA,YCRA,NCRA,IAAI,IAGI,NGPS)
        DIMENSION XCRA(NCRA),YCRA(NCRA),IAAI(NGPS),IAGI(NGPS)
        IAI1=-1
        DO 101 IGRP=1,NGPS
          IF (IAGI(IGRP).EQ.1) IAI1=IAAI(IGRP)
  101   CONTINUE
        IF (IAI1.GT.0) THEN
          IF (MAPACI(IAI1).EQ.1) THEN
            CALL GPL    (NCRA,XCRA,YCRA)
          END IF
        END IF
        RETURN
      END



      SUBROUTINE COLORB (XCRA,YCRA,NCRA,IAAI,IAGI,NGPS)
        DIMENSION XCRA(NCRA),YCRA(NCRA),IAAI(NGPS),IAGI(NGPS)
        IAI1=-1
        DO 101 IGRP=1,NGPS
          IF (IAGI(IGRP).EQ.1) IAI1=IAAI(IGRP)
  101   CONTINUE
        IF (IAI1.GT.0) THEN
          CALL GSFACI (MPISCI(IAI1)+100)
          CALL GFA    (NCRA-1,XCRA,YCRA)
        END IF
        RETURN
      END



      SUBROUTINE DRAWLB (XCRA,YCRA,NCRA,IAAI,IAGI,NGPS)
        DIMENSION XCRA(NCRA),YCRA(NCRA),IAAI(NGPS),IAGI(NGPS)
        IAI1=-1
        DO 101 IGRP=1,NGPS
          IF (IAGI(IGRP).EQ.1) IAI1=IAAI(IGRP)
  101   CONTINUE
        IF (IAI1.GT.0) THEN
          IF (MPIPAN(IAI1,'Water').NE.0) THEN
            CALL GPL    (NCRA,XCRA,YCRA)
          END IF
        END IF
        RETURN
      END



      SUBROUTINE MUPPER (CHRS)
        CHARACTER*(*) CHRS
        CHARACTER*26 LOWR,UPPR
        DATA LOWR / 'abcdefghijklmnopqrstuvwxyz' /
        DATA UPPR / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' /
        DO 102 I=1,LEN(CHRS)
          DO 101 J=1,26
            IF (CHRS(I:I).EQ.LOWR(J:J)) CHRS(I:I)=UPPR(J:J)
  101     CONTINUE
  102   CONTINUE
        RETURN
      END


      SUBROUTINE ZERODA (DPRA,LNTH)
        DOUBLE PRECISION DPRA(LNTH)
        DO 101 I=1,LNTH
          DPRA(I)=0.D0
  101   CONTINUE
        RETURN
      END


      SUBROUTINE DRSPGD (ISPG,IZN0,PLT0,PLN0,CLT0,CLN0,
     +                        IZN8,PLT8,PLN8,CLT8,CLN8,CHSZ)
C
        DIMENSION IZN0(139)
        DIMENSION PLT0(200,139),PLN0(200,139),CLT0(139),CLN0(139)
        DIMENSION IZN8(139)
        DIMENSION PLT8(200,139),PLN8(200,139),CLT8(139),CLN8(139)
C
C Declare a common block making it possible to set PLOTCHAR's internal
C parameter 'MA' = 4.
C
        COMMON /PCMP04/ PANG,PLAT,PLON
        SAVE   /PCMP04/
C
C Declare a character temporary.
C
        CHARACTER*5 CTMP
C
C Draw NAD1927 zone, if it's selected.
C
        IF (ISPG.EQ.0) THEN
C
          DO 103 ISPI=1,139
            IF (IZN0(ISPI).NE.0) THEN
              CALL SFLUSH
              CALL GSLWSC (6.)
              CALL GSPLCI (10)
              CALL MAPIT (PLT0(1,ISPI),PLN0(1,ISPI),0)
              DO 101 I=2,200
                CALL MAPIT (PLT0(I,ISPI),PLN0(I,ISPI),1)
  101         CONTINUE
              CALL MAPIT (PLT0(1,ISPI),PLN0(1,ISPI),2)
              CALL SFLUSH
              CALL GSLWSC (2.)
              CALL GSPLCI (12)
              CALL MAPIT (PLT0(1,ISPI),PLN0(1,ISPI),0)
              DO 102 I=2,200
                CALL MAPIT (PLT0(I,ISPI),PLN0(I,ISPI),1)
  102         CONTINUE
              CALL MAPIT (PLT0(1,ISPI),PLN0(1,ISPI),2)
            END IF
  103     CONTINUE
C
          CALL SFLUSH
          CALL GSLWSC (1.)
C
          CALL PCSETI ('MA',4)
          CALL PCSETR ('OR',1.E12)
C
          DO 104 ISPI=1,139
            IF (IZN0(ISPI).NE.0.AND.IZN0(ISPI).NE.1901) THEN
              WRITE (CTMP,'(I5)') 10000+IZN0(ISPI)
              PANG=0.
              PLAT=CLT0(ISPI)
              PLON=CLN0(ISPI)
              CALL PLCHHQ (0.,0.,CTMP(2:5),CHSZ,0.,0.)
            END IF
  104     CONTINUE
C
          CALL PCSETI ('MA',0)
          CALL PCSETR ('OR',0.)
C
          CALL SFLUSH
          CALL GSPLCI (1)
C
C Draw NAD1983 zone, if it's selected.
C
        ELSE IF (ISPG.EQ.8) THEN
C
          DO 107 ISPI=1,139
            IF (IZN8(ISPI).NE.0) THEN
              CALL SFLUSH
              CALL GSLWSC (6.)
              CALL GSPLCI (10)
              CALL MAPIT (PLT8(1,ISPI),PLN8(1,ISPI),0)
              DO 105 I=2,200
                CALL MAPIT (PLT8(I,ISPI),PLN8(I,ISPI),1)
  105         CONTINUE
              CALL MAPIT (PLT8(1,ISPI),PLN8(1,ISPI),2)
              CALL SFLUSH
              CALL GSLWSC (2.)
              CALL GSPLCI (12)
              CALL MAPIT (PLT8(1,ISPI),PLN8(1,ISPI),0)
              DO 106 I=2,200
                CALL MAPIT (PLT8(I,ISPI),PLN8(I,ISPI),1)
  106         CONTINUE
              CALL MAPIT (PLT8(1,ISPI),PLN8(1,ISPI),2)
            END IF
  107     CONTINUE
C
          CALL SFLUSH
          CALL GSLWSC (1.)
C
          CALL PCSETI ('MA',4)
          CALL PCSETR ('OR',1.E12)
C
          DO 108 ISPI=1,139
            IF (IZN8(ISPI).NE.0.AND.IZN8(ISPI).NE.1901) THEN
              WRITE (CTMP,'(I5)') 10000+IZN8(ISPI)
              PANG=0.
              PLAT=CLT8(ISPI)
              PLON=CLN8(ISPI)
              CALL PLCHHQ (0.,0.,CTMP(2:5),CHSZ,0.,0.)
            END IF
  108     CONTINUE
C
          CALL PCSETI ('MA',0)
          CALL PCSETR ('OR',0.)
C
          CALL SFLUSH
          CALL GSPLCI (1)
C
        END IF
C
        RETURN
C
      END


      SUBROUTINE INNPDP (IPDP,PTYP,PLAT,PLON,ROTA,LTYP,PLM1,PLM2,PLM3,
     +                   PLM4,IPRJ,IZON,ISPH,PARA,IDSL,ISPG,ODNM,ILVL,
     +                   CHSZ,LABL,ILLB)
C
        CHARACTER*(*) PTYP,LTYP,ODNM,LABL
        DIMENSION PLM1(2),PLM2(2),PLM3(2),PLM4(2)
        DOUBLE PRECISION PARA(15)
C
        PTYP='UT'
C
        PLAT=0.
        PLON=0.
        ROTA=0.
C
        LTYP='MA'
        PLM1(1)=0.
        PLM1(2)=0.
        PLM2(1)=0.
        PLM2(2)=0.
        PLM3(1)=0.
        PLM3(2)=0.
        PLM4(1)=0.
        PLM4(2)=0.
C
        IPRJ=0
        IZON=0
        ISPH=0
C
        IDSL=1
C
        ISPG=-1
C
        LABL=' '
C
        ILLB=0
C
        CALL ZERODA (PARA,15)
C
        PRINT * , ' '
C
        IF      (IPDP.EQ. 1) THEN
          PRINT * , 'Example 1 shows an original EZMAP transformation'
          PRINT * , 'called the Lambert Conformal Conic.'
          PTYP='LC'
          PLAT=55.
          ROTA=85.
        ELSE IF (IPDP.EQ. 2) THEN
          PRINT * , 'Example 2 shows an original EZMAP transformation'
          PRINT * , 'called the Gnomonic.'
          PTYP='GN'
          PLAT=40.
          PLON=-105.
        ELSE IF (IPDP.EQ. 3) THEN
          PRINT * , 'Example 3 shows an original EZMAP transformation'
          PRINT * , 'called the Mercator.'
          PTYP='ME'
          ILLB=-1
        ELSE IF (IPDP.EQ. 4) THEN
          PRINT * , 'Example 4 shows an original EZMAP transformation'
          PRINT * , 'called the Stereographic.'
          PTYP='ST'
        ELSE IF (IPDP.EQ. 5) THEN
          PRINT * , 'Example 5 shows an original EZMAP transformation'
          PRINT * , 'called the Lambert Equal-Area.'
          PTYP='LE'
        ELSE IF (IPDP.EQ. 6) THEN
          PRINT * , 'Example 6 shows an original EZMAP transformation'
          PRINT * , 'called the Orthographic.'
          PTYP='OR'
          PLAT=40.
          PLON=-105.
        ELSE IF (IPDP.EQ. 7) THEN
          PRINT * , 'Example 7 shows an original EZMAP transformation'
          PRINT * , 'called the Satellite View or Perspective.'
          PTYP='SV'
          PLAT=40.
          PLON=-105.
        ELSE IF (IPDP.EQ. 8) THEN
          PRINT * , 'Example 8 shows an original EZMAP transformation'
          PRINT * , 'called the Azimuthal Equidistant.'
          PTYP='AE'
        ELSE IF (IPDP.EQ. 9) THEN
          PRINT * , 'Example 9 shows an original EZMAP transformation'
          PRINT * , 'called the Robinson, which was developed by the'
          PRINT * , 'National Geographic Society for use in world maps.'
          PTYP='RO'
        ELSE IF (IPDP.EQ.10) THEN
          PRINT * , 'Example 10 shows an original EZMAP transformation'
          PRINT * , 'called the Mollweide.'
          PTYP='MO'
        ELSE IF (IPDP.EQ.11) THEN
          PRINT * , 'Example 11 shows an original EZMAP transformation'
          PRINT * , 'called the Cylindrical Equidistant.'
          PTYP='CE'
          ILLB=-1
        ELSE IF (IPDP.EQ.12) THEN
          PRINT * , 'Example 12 shows an original EZMAP transformation'
          PRINT * , 'called the Cylindrical Equal-Area.'
          PTYP='EA'
          ILLB=-1
        ELSE IF (IPDP.EQ.13) THEN
          PRINT * , 'Example 13 shows an original EZMAP transformation'
          PRINT * , 'called the Rotated Mercator.'
          PTYP='RM'
          ROTA=45.
          LTYP='CO'
          PLM1(1)=   2.
          PLM2(1)=-135.
          PLM3(1)=  37.
          PLM4(1)=  67.
        ELSE IF (IPDP.EQ.14) THEN
          PRINT * , 'Example 14 shows a USGS transformation called the'
          PRINT * , 'UTM system (UTM="Universal Transverse Mercator"),'
          PRINT * , 'in the Northern Hemisphere.  The area shown is'
          PRINT * , 'that for which distortion is minimal, on either'
          PRINT * , 'side of a particular meridian; normally, the UTM'
          PRINT * , 'system is used for large-scale maps which are more'
          PRINT * , 'nearly square.'
          IPRJ=1
          IZON=13
          IDSL=4
        ELSE IF (IPDP.EQ.15) THEN
          PRINT * , 'Example 15 shows a USGS transformation called the'
          PRINT * , 'UTM system (UTM="Universal Transverse Mercator"),'
          PRINT * , 'in the Southern Hemisphere.  The area shown is'
          PRINT * , 'that for which distortion is minimal, on either'
          PRINT * , 'side of a particular meridian; normally, the UTM'
          PRINT * , 'system is used for large-scale maps which are more'
          PRINT * , 'nearly square.'
          IPRJ=1
          IZON=-35
          IDSL=4
        ELSE IF (IPDP.EQ.16) THEN
          PRINT * , 'Example 16 shows a USGS transformation called the'
          PRINT * , 'State Plane coordinate system.  The zone shown is'
          PRINT * , 'that for New Hampshire.'
          IPRJ=2
          IZON=2800
          ISPH=0
          IDSL=4
          ILLB=-1
        ELSE IF (IPDP.EQ.17) THEN
          PRINT * , 'Example 17 shows a USGS transformation called the'
          PRINT * , 'State Plane coordinate system.  The zone shown is'
          PRINT * , 'that for Nebraska.'
          IPRJ=2
          IZON=2600
          ISPH=8
          IDSL=4
          ILLB=-1
        ELSE IF (IPDP.EQ.18) THEN
          PRINT * , 'Example 18 shows a USGS transformation called the'
          PRINT * , 'State Plane coordinate system.  The zone shown is'
          PRINT * , 'that for the Alaskan panhandle.'
          IPRJ=2
          IZON=5001
          ISPH=8
          IDSL=4
          ILLB=-1
        ELSE IF (IPDP.EQ.19) THEN
          PRINT * , 'Example 19 shows a USGS transformation called the'
          PRINT * , 'Alber''s Equal-Area Conic, in the Northern'
          PRINT * , 'Hemisphere.'
          IPRJ=3
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=         50.D0  !  LATITUDE OF 1ST STANDARD PARALLEL
          PARA( 4)=         60.D0  !  LATITUDE OF 2ND STANDARD PARALLEL
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=         90.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.20) THEN
          PRINT * , 'Example 20 shows a USGS transformation called the'
          PRINT * , 'Alber''s Equal-Area Conic, in the Southern'
          PRINT * , 'Hemisphere.'
          IPRJ=3
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=        -50.D0  !  LATITUDE OF 1ST STANDARD PARALLEL
          PARA( 4)=        -60.D0  !  LATITUDE OF 2ND STANDARD PARALLEL
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=        -90.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.21) THEN
          PRINT * , 'Example 21 shows a USGS transformation called the'
          PRINT * , 'Lambert Conformal Conic, with two standard'
          PRINT * , 'parallels, in the Northern Hemisphere.'
          IPRJ=4
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=         50.D0  !  LATITUDE OF 1ST STANDARD PARALLEL
          PARA( 4)=         60.D0  !  LATITUDE OF 2ND STANDARD PARALLEL
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=         55.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.22) THEN
          PRINT * , 'Example 22 shows a USGS transformation called the'
          PRINT * , 'Lambert Conformal Conic, with two standard'
          PRINT * , 'parallels, in the Southern Hemisphere.'
          IPRJ=4
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=        -50.D0  !  LATITUDE OF 1ST STANDARD PARALLEL
          PARA( 4)=        -60.D0  !  LATITUDE OF 2ND STANDARD PARALLEL
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=        -55.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.23) THEN
          PRINT * , 'Example 23 shows a USGS transformation called the'
          PRINT * , 'Mercator.'
          IPRJ=5
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=          0.D0  !  LATITUDE OF TRUE SCALE
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.24) THEN
          PRINT * , 'Example 24 shows a USGS transformation called the'
          PRINT * , 'Polar Stereographic, in the Northern Hemisphere.'
          IPRJ=6
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 5)=       -105.D0  !  LONGITUDE DOWN FROM NORTH, UP FROM SOUTH
          PARA( 6)=        +45.D0  !  LATITUDE OF TRUE SCALE
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.25) THEN
          PRINT * , 'Example 25 shows a USGS transformation called the'
          PRINT * , 'Polar Stereographic, in the Southern Hemisphere.'
          IPRJ=6
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 5)=       -105.D0  !  LONGITUDE DOWN FROM NORTH, UP FROM SOUTH
          PARA( 6)=        -45.D0  !  LATITUDE OF TRUE SCALE
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.26) THEN
          PRINT * , 'Example 26 shows a USGS transformation called the'
          PRINT * , 'Polyconic.'
          IPRJ=7
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=          0.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.27) THEN
          PRINT * , 'Example 27 shows a USGS transformation called the'
          PRINT * , 'Equidistant Conic, with one standard parallel, in'
          PRINT * , 'the Northern Hemisphere.'
          IPRJ=8
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=         40.D0  !  LATITUDE OF STANDARD PARALLEL
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=         90.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
          PARA( 9)=          0.D0  !  ZERO
        ELSE IF (IPDP.EQ.28) THEN
          PRINT * , 'Example 28 shows a USGS transformation called the'
          PRINT * , 'Equidistant Conic, with one standard parallel, in'
          PRINT * , 'the Southern Hemisphere.'
          IPRJ=8
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=        -40.D0  !  LATITUDE OF STANDARD PARALLEL
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=        -90.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
          PARA( 9)=          0.D0  !  ZERO
        ELSE IF (IPDP.EQ.29) THEN
          PRINT * , 'Example 29 shows a USGS transformation called the'
          PRINT * , 'Equidistant Conic, with two standard parallels, in'
          PRINT * , 'the Northern Hemisphere.'
          IPRJ=8
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=         40.D0  !  LATITUDE OF 1ST STANDARD PARALLEL
          PARA( 4)=         55.D0  !  LATITUDE OF 2ND STANDARD PARALLEL
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=         90.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
          PARA( 9)=          1.D0  !  NON-ZERO
        ELSE IF (IPDP.EQ.30) THEN
          PRINT * , 'Example 30 shows a USGS transformation called the'
          PRINT * , 'Equidistant Conic, with two standard parallels, in'
          PRINT * , 'the Southern Hemisphere.'
          IPRJ=8
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=        -40.D0  !  LATITUDE OF 1ST STANDARD PARALLEL
          PARA( 4)=        -55.D0  !  LATITUDE OF 2ND STANDARD PARALLEL
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=        -90.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
          PARA( 9)=          1.D0  !  NON-ZERO
        ELSE IF (IPDP.EQ.31) THEN
          PRINT * , 'Example 31 shows a USGS transformation called the'
          PRINT * , 'Transverse Mercator.  The area shown is that for'
          PRINT * , 'which distortion is minimal, on either side of a'
          PRINT * , 'particular meridian; normally, this system is used'
          PRINT * , 'for large-scale maps which are more nearly square.'
          IPRJ=9
          ISPH=-1
          PARA( 1)=   6378206.4D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=   6356583.8D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=          1.D0  !  SCALE FACTOR AT CENTRAL MERIDIAN
          PARA( 5)=         25.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=          0.D0  !  LATITUDE OF ORIGIN
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.32) THEN
          PRINT * , 'Example 32 shows a USGS transformation called the'
          PRINT * , 'Stereographic.'
          IPRJ=10
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTER OF PROJECTION
          PARA( 6)=         40.D0  !  LATITUDE OF CENTER OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.33) THEN
          PRINT * , 'Example 33 shows a USGS transformation called the'
          PRINT * , 'Lambert Azimuthal Equal-Area.'
          IPRJ=11
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTER OF PROJECTION
          PARA( 6)=         40.D0  !  LATITUDE OF CENTER OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.34) THEN
          PRINT * , 'Example 34 shows a USGS transformation: the'
          PRINT * , 'Azimuthal Equidistant.'
          IPRJ=12
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTER OF PROJECTION
          PARA( 6)=         40.D0  !  LATITUDE OF CENTER OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.35) THEN
          PRINT * , 'Example 35 shows a USGS transformation called the'
          PRINT * , 'Gnomonic.'
          IPRJ=13
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTER OF PROJECTION
          PARA( 6)=         40.D0  !  LATITUDE OF CENTER OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.36) THEN
          PRINT * , 'Example 36 shows a USGS transformation called the'
          PRINT * , 'Orthographic.'
          IPRJ=14
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTER OF PROJECTION
          PARA( 6)=         40.D0  !  LATITUDE OF CENTER OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.37) THEN
          PRINT * , 'Example 37 shows a USGS transformation called the'
          PRINT * , 'Perspective (the same as the Satellite-View, but'
          PRINT * , 'only allows for a vertical view).'
          IPRJ=15
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 3)=          3.D0  !  HEIGHT OF PERSPECTIVE POINT ABOVE SURFACE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTER OF PROJECTION
          PARA( 6)=         40.D0  !  LATITUDE OF CENTER OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.38) THEN
          PRINT * , 'Example 38 shows a USGS transformation called the'
          PRINT * , 'Sinusoidal.'
          IPRJ=16
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.39) THEN
          PRINT * , 'Example 39 shows a USGS transformation called the'
          PRINT * , 'Equirectangular (same as Cylindrical Equidistant).'
          IPRJ=17
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 6)=          0.D0  !  LATITUDE OF TRUE SCALE
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.40) THEN
          PRINT * , 'Example 40 shows a USGS transformation called the'
          PRINT * , 'Miller Cylindrical.'
          IPRJ=18
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.41) THEN
          PRINT * , 'Example 41 shows a USGS transformation called the'
          PRINT * , 'Van der Grinten I.'
          IPRJ=19
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.42) THEN
          PRINT * , 'Example 42 shows a USGS transformation called the'
          PRINT * , 'Oblique Mercator.  The projection''s center line'
          PRINT * , 'is specified by giving two points that project'
          PRINT * , 'into it.  Currently, there is a bug which prevents'
          PRINT * , 'this specification method from working properly.'
          IPRJ=20
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=          1.D0  !  SCALE FACTOR AT CENTER OF PROJECTION
          PARA( 6)=         40.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
          PARA( 9)=       -100.D0  !  LONGITUDE, POINT 1
          PARA(10)=         30.D0  !  LATITUDE, POINT 1
          PARA(11)=        -80.D0  !  LONGITUDE, POINT 2
          PARA(12)=         50.D0  !  LATITUDE, POINT 2
          PARA(13)=          0.D0  !  ZERO
        ELSE IF (IPDP.EQ.43) THEN
          PRINT * , 'Example 43 shows a USGS transformation called the'
          PRINT * , 'Oblique Mercator.  The projection''s center line'
          PRINT * , 'is specified by specifying its azimuth angle.'
          IPRJ=20
          ISPH=-1
          PARA( 1)=          1.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=        .007D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=          1.D0  !  SCALE FACTOR AT CENTER OF PROJECTION
          PARA( 4)=         45.D0  !  AZIMUTH ANGLE OF CENTER LINE (E OF N)
          PARA( 5)=       -105.D0  !  LONGITUDE OF ORIGIN OF PROJECTION
          PARA( 6)=         40.D0  !  LATITUDE OF ORIGIN OF PROJECTION
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
          PARA(13)=          1.D0  !  NON-ZERO
        ELSE IF (IPDP.EQ.44) THEN
          PRINT * , 'Example 44 shows a USGS transformation called the'
          PRINT * , 'Robinson, which was developed by the National'
          PRINT * , 'Geographic Society for use in world maps.'
          IPRJ=21
          ISPH=-1
          PARA( 1)=          1.D0  !  RADIUS OF SPHERE OF REFERENCE
          PARA( 5)=       -105.D0  !  LONGITUDE OF CENTRAL MERIDIAN
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.45) THEN
          PRINT * , 'Example 45 shows a USGS transformation called the'
          PRINT * , 'Space Oblique Mercator.  Each of the views that'
          PRINT * , 'one can select shows the area passed over by a'
          PRINT * , 'chosen Landsat satellite during a particular'
          PRINT * , 'orbit.'
          IPRJ=22
          PARA( 1)=          0.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=          0.D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 3)=          1.D0  !  LANDSAT NUMBER (1-5)
          PARA( 4)=         84.D0  !  PATH (1-251 for LN 1-3, 1-233 for LN 4-5)
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
        ELSE IF (IPDP.EQ.46) THEN
          PRINT * , 'Example 46 shows a USGS transformation called the'
          PRINT * , 'Modified Stereographic for Alaska.  It is centered'
          PRINT * , 'at Lat 64N, Long 152W and stretched in such a way'
          PRINT * , 'as to allow all of Alaska to be shown on a single'
          PRINT * , 'map with relatively little distortion.'
          IPRJ=23
          PARA( 1)=          0.D0  !  SEMIMAJOR AXIS OF ELLIPSOID
          PARA( 2)=          0.D0  !  SQUARE OF ELLIPSOID'S ECCENTRICITY
          PARA( 7)=          0.D0  !  FALSE EASTING
          PARA( 8)=          0.D0  !  FALSE NORTHING
          IDSL=4
          ILLB=-1
        ELSE IF (IPDP.EQ.47) THEN
          PRINT * , 'Example 47 shows the State Plane zones for the'
          PRINT * , 'northeastern states of the lower 48, using the'
          PRINT * , 'North American Datum of 1927.'
          PTYP='ST'
          PLAT=41.
          PLON=-76.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=0
          CHSZ=.24
          LABL='STATE PLANE ZONES - NORTHEAST - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.48) THEN
          PRINT * , 'Example 48 shows the State Plane zones for the'
          PRINT * , 'northeastern states of the lower 48, using the'
          PRINT * , 'North American Datum of 1983.'
          PTYP='ST'
          PLAT=41.
          PLON=-76.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=8
          CHSZ=.24
          LABL='STATE PLANE ZONES - NORTHEAST - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.49) THEN
          PRINT * , 'Example 49 shows the State Plane zones for the'
          PRINT * , 'southeastern states of the lower 48, using the'
          PRINT * , 'North American Datum of 1927.'
          PTYP='ST'
          PLAT=31.
          PLON=-84.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=0
          CHSZ=.24
          LABL='STATE PLANE ZONES - SOUTHEAST - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.50) THEN
          PRINT * , 'Example 50 shows the State Plane zones for the'
          PRINT * , 'southeastern states of the lower 48, using the'
          PRINT * , 'North American Datum of 1983.'
          PTYP='ST'
          PLAT=31.
          PLON=-84.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=8
          CHSZ=.24
          LABL='STATE PLANE ZONES - SOUTHEAST - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.51) THEN
          PRINT * , 'Example 51 shows the State Plane zones for the'
          PRINT * , 'northeast central states of the lower 48, using'
          PRINT * , 'the North American Datum of 1927.'
          PTYP='ST'
          PLAT=42.
          PLON=-92.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=0
          CHSZ=.24
          LABL='STATE PLANE ZONES - NORTHEAST CENTRAL - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.52) THEN
          PRINT * , 'Example 52 shows the State Plane zones for the'
          PRINT * , 'northeast central states of the lower 48, using'
          PRINT * , 'the North American Datum of 1983.'
          PTYP='ST'
          PLAT=42.
          PLON=-92.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=8
          CHSZ=.24
          LABL='STATE PLANE ZONES - NORTHEAST CENTRAL - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.53) THEN
          PRINT * , 'Example 53 shows the State Plane zones for the'
          PRINT * , 'southeast central states of the lower 48, using'
          PRINT * , 'the North American Datum of 1927.'
          PTYP='ST'
          PLAT=33.
          PLON=-92.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=0
          CHSZ=.24
          LABL='STATE PLANE ZONES - SOUTHEAST CENTRAL - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.54) THEN
          PRINT * , 'Example 54 shows the State Plane zones for the'
          PRINT * , 'southeast central states of the lower 48, using'
          PRINT * , 'the North American Datum of 1983.'
          PTYP='ST'
          PLAT=33.
          PLON=-92.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=8
          CHSZ=.24
          LABL='STATE PLANE ZONES - SOUTHEAST CENTRAL - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.55) THEN
          PRINT * , 'Example 55 shows the State Plane zones for the'
          PRINT * , 'northwest central states of the lower 48, using'
          PRINT * , 'the North American Datum of 1927.'
          PTYP='ST'
          PLAT=42.
          PLON=-104.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=0
          CHSZ=.24
          LABL='STATE PLANE ZONES - NORTHWEST CENTRAL - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.56) THEN
          PRINT * , 'Example 56 shows the State Plane zones for the'
          PRINT * , 'northwest central states of the lower 48, using'
          PRINT * , 'the North American Datum of 1983.'
          PTYP='ST'
          PLAT=42.
          PLON=-104.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=8
          CHSZ=.24
          LABL='STATE PLANE ZONES - NORTHWEST CENTRAL - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.57) THEN
          PRINT * , 'Example 57 shows the State Plane zones for the'
          PRINT * , 'southwest central states of the lower 48, using'
          PRINT * , 'the North American Datum of 1927.'
          PTYP='ST'
          PLAT=33.
          PLON=-104.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=0
          CHSZ=.24
          LABL='STATE PLANE ZONES - SOUTHWEST CENTRAL - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.58) THEN
          PRINT * , 'Example 58 shows the State Plane zones for the'
          PRINT * , 'southwest central states of the lower 48, using'
          PRINT * , 'the North American Datum of 1983.'
          PTYP='ST'
          PLAT=33.
          PLON=-104.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=8
          CHSZ=.24
          LABL='STATE PLANE ZONES - SOUTHWEST CENTRAL - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.59) THEN
          PRINT * , 'Example 59 shows the State Plane zones for the'
          PRINT * , 'northwestern states of the lower 48, using the'
          PRINT * , 'North American Datum of 1927.'
          PTYP='ST'
          PLAT=42.
          PLON=-116.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=0
          CHSZ=.24
          LABL='STATE PLANE ZONES - NORTHWEST - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.60) THEN
          PRINT * , 'Example 60 shows the State Plane zones for the'
          PRINT * , 'northwestern states of the lower 48, using the'
          PRINT * , 'North American Datum of 1983.'
          PTYP='ST'
          PLAT=42.
          PLON=-116.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=8
          CHSZ=.24
          LABL='STATE PLANE ZONES - NORTHWEST - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.61) THEN
          PRINT * , 'Example 61 shows the State Plane zones for the'
          PRINT * , 'southwestern states of the lower 48, using the'
          PRINT * , 'North American Datum of 1927.'
          PTYP='ST'
          PLAT=38.
          PLON=-116.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=0
          CHSZ=.24
          LABL='STATE PLANE ZONES - SOUTHWEST - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.62) THEN
          PRINT * , 'Example 62 shows the State Plane zones for the'
          PRINT * , 'southwestern states of the lower 48, using the'
          PRINT * , 'North American Datum of 1983.'
          PTYP='ST'
          PLAT=38.
          PLON=-116.
          LTYP='AN'
          PLM1(1)=8.
          PLM2(1)=8.
          PLM3(1)=8.
          PLM4(1)=8.
          IDSL=4
          ISPG=8
          CHSZ=.24
          LABL='STATE PLANE ZONES - SOUTHWEST - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.63) THEN
          PRINT * , 'Example 63 shows the State Plane zones for Alaska,'
          PRINT * , 'using the North American Datum of 1927.'
          PTYP='ST'
          PLAT=64.
          PLON=-152.
          LTYP='AN'
          PLM1(1)=23.
          PLM2(1)=15.
          PLM3(1)=19.
          PLM4(1)=19.
          IDSL=4
          ISPG=0
          CHSZ=.75
          LABL='STATE PLANE ZONES - ALASKA - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.64) THEN
          PRINT * , 'Example 64 shows the State Plane zones for Alaska,'
          PRINT * , 'using the North American Datum of 1983.'
          PTYP='ST'
          PLAT=64.
          PLON=-152.
          LTYP='AN'
          PLM1(1)=23.
          PLM2(1)=15.
          PLM3(1)=19.
          PLM4(1)=19.
          IDSL=4
          ISPG=8
          CHSZ=.75
          LABL='STATE PLANE ZONES - ALASKA - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.65) THEN
          PRINT * , 'Example 65 shows the State Plane zones for Hawaii,'
          PRINT * , 'using the North American Datum of 1927.'
          PTYP='ST'
          PLAT=21.
          PLON=-157.5
          LTYP='AN'
          PLM1(1)=3.
          PLM2(1)=3.
          PLM3(1)=3.
          PLM4(1)=3.
          IDSL=4
          ISPG=0
          CHSZ=.075
          LABL='STATE PLANE ZONES - HAWAII - NAD1927'
          ILLB=-1
        ELSE IF (IPDP.EQ.66) THEN
          PRINT * , 'Example 66 shows the State Plane zones for Hawaii,'
          PRINT * , 'using the North American Datum of 1927.'
          PTYP='ST'
          PLAT=21.
          PLON=-157.5
          LTYP='AN'
          PLM1(1)=3.
          PLM2(1)=3.
          PLM3(1)=3.
          PLM4(1)=3.
          IDSL=4
          ISPG=8
          CHSZ=.075
          LABL='STATE PLANE ZONES - HAWAII - NAD1983'
          ILLB=-1
        ELSE IF (IPDP.EQ.67) THEN
          PRINT * , 'Example 67 shows all the State Plane zones, using'
          PRINT * , 'the North American Datum of 1927.'
          PTYP='LE'
          PLAT=40.
          PLON=-140.
          LTYP='AN'
          PLM1(1)=70.
          PLM2(1)=70.
          PLM3(1)=70.
          PLM4(1)=70.
          IDSL=4
          ISPG=0
          CHSZ=1.75
          LABL='STATE PLANE ZONES - ALL - NAD1927'
        ELSE IF (IPDP.EQ.68) THEN
          PRINT * , 'Example 68 shows all the State Plane zones, using'
          PRINT * , 'the North American Datum of 1983.'
          PTYP='LE'
          PLAT=45.
          PLON=-112.
          LTYP='AN'
          PLM1(1)=48.
          PLM2(1)=48.
          PLM3(1)=48.
          PLM4(1)=48.
          IDSL=4
          ISPG=8
          CHSZ=1.7
          LABL='STATE PLANE ZONES - ALL - NAD1983'
        END IF
C
        IF (IPDP.GE.47.AND.IPDP.LE.68) THEN
          PRINT * , ' '
          PRINT * , 'The area for which each of the State Plane zones'
          PRINT * , 'shown is meant to be used is shown as a rectangle'
          PRINT * , 'with rounded corners, so that one can more easily'
          PRINT * , 'sort out how the zones overlap.  Each zone number'
          PRINT * , 'is approximately centered within its zone.'
        END IF
C
        IF (IDSL.EQ.1) THEN
          IF (ODNM.NE.'E1'.AND.ODNM.NE.'E2'.AND.ODNM.NE.'E3') THEN
            ODNM='CO'
          ELSE
            ILVL=1
          END IF
        ELSE
          IF (ODNM.NE.'E1'.AND.ODNM.NE.'E2'.AND.ODNM.NE.'E3') THEN
            ODNM='PS'
          ELSE
            ILVL=4
          END IF
        END IF
C
        IPDP=MOD(IPDP,68)+1
C
        RETURN
C
      END


      SUBROUTINE DRPTLB (PTYP,IPRJ,IZON,ISPH,PARA,IPRF,ODNM,ILVL)
C
        CHARACTER*(*) PTYP,ODNM
        DOUBLE PRECISION PARA(15)
C
        CHARACTER*82 PTLB
C
        PTLB=' '
C
        IF      (PTYP.EQ.'LC') THEN
          PTLB='LAMBERT CONFORMAL CONIC'
        ELSE IF (PTYP.EQ.'GN') THEN
          PTLB='GNOMONIC'
        ELSE IF (PTYP.EQ.'ME') THEN
          PTLB='MERCATOR'
        ELSE IF (PTYP.EQ.'RM') THEN
          PTLB='ROTATED MERCATOR'
        ELSE IF (PTYP.EQ.'ST') THEN
          PTLB='STEREOGRAPHIC'
        ELSE IF (PTYP.EQ.'LE') THEN
          PTLB='LAMBERT EQUAL-AREA'
        ELSE IF (PTYP.EQ.'OR') THEN
          PTLB='ORTHOGRAPHIC'
        ELSE IF (PTYP.EQ.'SV') THEN
          PTLB='SATELLITE-VIEW'
        ELSE IF (PTYP.EQ.'AE') THEN
          PTLB='AZIMUTHAL EQUIDISTANT'
        ELSE IF (PTYP.EQ.'RO') THEN
          PTLB='ROBINSON'
        ELSE IF (PTYP.EQ.'MO') THEN
          PTLB='MOLLWEIDE'
        ELSE IF (PTYP.EQ.'CE') THEN
          PTLB='CYLINDRICAL EQUIDISTANT'
        ELSE IF (PTYP.EQ.'EA') THEN
          PTLB='CYLINDRICAL EQUAL-AREA'
        ELSE IF (PTYP.EQ.'UT') THEN
          IF (IPRJ.EQ.1) THEN
            IF (IZON.GT.0) THEN
              PTLB='USGS - UTM COORDINATES - Zone 00 - North'
            ELSE
              PTLB='USGS - UTM COORDINATES - Zone 00 - South'
            END IF
            CALL PUTINT (PTLB(31:32),IZON)
          ELSE IF (IPRJ.EQ.2) THEN
            PTLB( 1:38)='USGS - STATE PLANE COORDINATES (XX) - '
            PTLB(39:57)='NAD0000 - Zone 0000'
            IF (IPRF.EQ.4) THEN
              PTLB(33:34)='LC'
            ELSE IF (IPRF.EQ.7) THEN
              PTLB(33:34)='PC'
            ELSE IF (IPRF.EQ.9) THEN
              PTLB(33:34)='TM'
            ELSE IF (IPRF.EQ.20) THEN
              PTLB(33:34)='OM'
            END IF
            IF (ISPH.EQ.0) THEN
              PTLB(42:45)='1927'
            ELSE
              PTLB(42:45)='1983'
            END IF
            CALL PUTINT (PTLB(54:57),IZON)
          ELSE IF (IPRJ.EQ.3) THEN
            IF (PARA(3)+PARA(4).GT.0.D0) THEN
              PTLB='USGS - ALBER''S EQUAL-AREA CONIC - North'
            ELSE
              PTLB='USGS - ALBER''S EQUAL-AREA CONIC - South'
            END IF
          ELSE IF (IPRJ.EQ.4) THEN
            IF (PARA(3)+PARA(4).GT.0.D0) THEN
                PTLB='USGS - LAMBERT CONFORMAL CONIC - North'
            ELSE
                PTLB='USGS - LAMBERT CONFORMAL CONIC - South'
            END IF
          ELSE IF (IPRJ.EQ.5) THEN
            PTLB='USGS - MERCATOR'
          ELSE IF (IPRJ.EQ.6) THEN
            IF (PARA(6).GT.0.D0) THEN
              PTLB='USGS - POLAR STEREOGRAPHIC - North'
            ELSE
              PTLB='USGS - POLAR STEREOGRAPHIC - South'
            END IF
          ELSE IF (IPRJ.EQ.7) THEN
            PTLB='USGS - POLYCONIC'
          ELSE IF (IPRJ.EQ.8) THEN
            IF (PARA(9).EQ.0.D0) THEN
              IF (PARA(3).GT.0.D0) THEN
                PTLB='USGS - EQUIDISTANT CONIC - Type A - North'
              ELSE
                PTLB='USGS - EQUIDISTANT CONIC - Type A - South'
              END IF
            ELSE
              IF (PARA(3)+PARA(4).GT.0.D0) THEN
                PTLB='USGS - EQUIDISTANT CONIC - Type B - North'
              ELSE
                PTLB='USGS - EQUIDISTANT CONIC - Type B - South'
              END IF
            END IF
          ELSE IF (IPRJ.EQ.9) THEN
            PTLB='USGS - TRANSVERSE MERCATOR'
          ELSE IF (IPRJ.EQ.10) THEN
            PTLB='USGS - STEREOGRAPHIC'
          ELSE IF (IPRJ.EQ.11) THEN
            PTLB='USGS - LAMBERT EQUAL-AREA'
          ELSE IF (IPRJ.EQ.12) THEN
            PTLB='USGS - AZIMUTHAL EQUIDISTANT'
          ELSE IF (IPRJ.EQ.13) THEN
            PTLB='USGS - GNOMONIC'
          ELSE IF (IPRJ.EQ.14) THEN
            PTLB='USGS - ORTHOGRAPHIC'
          ELSE IF (IPRJ.EQ.15) THEN
            PTLB='USGS - PERSPECTIVE'
          ELSE IF (IPRJ.EQ.16) THEN
            PTLB='USGS - SINUSOIDAL'
          ELSE IF (IPRJ.EQ.17) THEN
            PTLB='USGS - EQUIRECTANGULAR'
          ELSE IF (IPRJ.EQ.18) THEN
            PTLB='USGS - MILLER CYLINDRICAL'
          ELSE IF (IPRJ.EQ.19) THEN
            PTLB='USGS - VAN DER GRINTEN I'
          ELSE IF (IPRJ.EQ.20) THEN
            IF (PARA(13).EQ.0.D0) THEN
              PTLB='USGS - OBLIQUE MERCATOR - Type A'
            ELSE
              PTLB='USGS - OBLIQUE MERCATOR - Type B'
            END IF
          ELSE IF (IPRJ.EQ.21) THEN
            PTLB='USGS - ROBINSON'
          ELSE IF (IPRJ.EQ.22) THEN
            PTLB='USGS - SPACE OBLIQUE MERCATOR'
          ELSE IF (IPRJ.EQ.23) THEN
            PTLB='USGS - MODIFIED STEREOGRAPHIC FOR ALASKA'
          END IF
        END IF
C
C Add the dataset specifier.
C
        ILST=MPILNB(PTLB)
        PTLB(ILST+ 1:ILST+12)=' - Database '
        PTLB(ILST+13:ILST+14)=ODNM(1:2)
        ILST=ILST+14
        IF (ODNM.EQ.'E1'.OR.ODNM.EQ.'E2'.OR.ODNM.EQ.'E3') THEN
          PTLB(ILST+1:ILST+10)=' at level '
          CALL PUTINT (PTLB(ILST+11:ILST+11),ILVL)
          ILST=ILST+11
        END IF
C
        IF (PTLB.NE.' ') CALL PLCHHQ (CFUX(.5),CFUY(.05),PTLB(1:ILST),
     +                                                     .012,0.,0.)
C
        RETURN
C
      END


      SUBROUTINE PUTINT (CHRS,IVAL)
        CHARACTER*(*) CHRS
        CHARACTER*10 DGTS
        DATA DGTS / '0123456789' /
        IIII=ABS(IVAL)
        DO 101 I=LEN(CHRS),1,-1
          IDIG=MOD(IIII,10)+1
          CHRS(I:I)=DGTS(IDIG:IDIG)
          IIII=IIII/10
  101   CONTINUE
        RETURN
      END


      SUBROUTINE DRAWSP (ICLR,PLAT,PLON,NOTA)
C
C Draw a spiral, in the color specified by ICLR, which starts at the
C point with latitude PLAT and longitude PLON and ends at the point
C opposite it on the globe.  NOTA is the number of times the spiral
C winds around the globe.
C
        NPTS=720*NOTA-1
C
        RDIV=REAL(NPTS+1)
C
        CALL SFLUSH
        CALL GSPLCI (ICLR)
C
        DO 101 IPNT=0,NPTS+1
          RLAT=-90.+180.*REAL(IPNT)/RDIV
          RLON=.5*REAL(MOD(IPNT,720))
          CALL ROTATE (2,-PLAT-90.,RLAT,RLON)
          CALL ROTATE (3, PLON    ,RLAT,RLON)
          IF (IPNT.EQ.0) THEN
            CALL MAPIT (RLAT,RLON,0)
          ELSE IF (IPNT.LE.NPTS) THEN
            CALL MAPIT (RLAT,RLON,1)
          ELSE
            CALL MAPIT (RLAT,RLON,2)
          END IF
  101   CONTINUE
C
        CALL SFLUSH
        CALL GSPLCI (1)
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE ROTATE (IAXS,ANGL,RLAT,RLON)
C
C Define multiplicative constants to convert from degrees to radians
C and from radians to degrees.
C
        DATA DTOR / .017453292519943 /
        DATA RTOD / 57.2957795130823 /
C
C This routine rotates the point with latitude RLAT and longitude RLON
C by the angle ANGL about the axis specified by IAXS (1 for the U axis,
C 2 for the V axis, 3 for the W axis).  A right-handed coordinate system
C is assumed.
C
        UCRD=COS(DTOR*RLAT)*COS(DTOR*RLON)
        VCRD=COS(DTOR*RLAT)*SIN(DTOR*RLON)
        WCRD=SIN(DTOR*RLAT)
C
        SINA=SIN(DTOR*ANGL)
        COSA=COS(DTOR*ANGL)
C
        UTMP=UCRD
        VTMP=VCRD
        WTMP=WCRD
C
        IF (IAXS.EQ.1) THEN
          VCRD=VTMP*COSA-WTMP*SINA
          WCRD=WTMP*COSA+VTMP*SINA
        ELSE IF (IAXS.EQ.2) THEN
          UCRD=UTMP*COSA+WTMP*SINA
          WCRD=WTMP*COSA-UTMP*SINA
        ELSE
          UCRD=UTMP*COSA-VTMP*SINA
          VCRD=VTMP*COSA+UTMP*SINA
        END IF
C
        RLAT=RTOD*ASIN(WCRD)
        RLON=RTOD*ATAN2(VCRD,UCRD)
C
        RETURN
C
      END


      SUBROUTINE CONINV (ICLR,NPTS)
C
C Declare one of the USGS code's common blocks to gain access to its
C error flag.
C
        COMMON /ERRMZ0/ IERR
        SAVE   /ERRMZ0/
C
C Declare the arrays required for contouring values returned by MAPTRI.
C
        PARAMETER (MPTS=501,LRWK=10000,LIWK=10000)
C
        DIMENSION RLAT(MPTS,MPTS),RLON(MPTS,MPTS),RWRK(LRWK),IWRK(LIWK)
C
        CALL CPSETR ('VPL - VIEWPORT LEFT  ',0.)
        CALL CPSETR ('VPR - VIEWPORT RIGHT ',1.)
        CALL CPSETR ('VPB - VIEWPORT BOTTOM',0.)
        CALL CPSETR ('VPT - VIEWPORT TOP   ',1.)
        CALL CPSETI ('CLS - CONTOUR LEVEL SELECTOR',1)
        CALL CPSETI ('CIS - CONTOUR INTERVAL SELECTOR',5)
        CALL CPSETI ('LIS - LABEL INTERVAL SELECTOR',100)
        CALL CPSETR ('SPV - SPECIAL VALUE',1.E12)
        CALL CPSETI ('RWC - REAL WORKSPACE FOR CONTOURS',5000)
        CALL CPSETI ('PAI - PARAMETER ARRAY INDEX',-2)
        CALL CPSETI ('CLU - CONTOUR LEVEL USE FLAG',1)
C
        CALL SFLUSH
        CALL GSPLCI (ICLR)
        CALL GSTXCI (ICLR)
        CALL GSMKSC (50./REAL(NPTS-1))
C
        DO 102 I=1,NPTS
          XVAL=CFUX(REAL(I-1)/REAL(NPTS-1))
          DO 101 J=1,NPTS
            YVAL=CFUY(REAL(J-1)/REAL(NPTS-1))
            IERR=0
            CALL MAPTRI (XVAL,YVAL,RLAT(I,J),RLON(I,J))
            IF (RLAT(I,J).EQ.1.E12) THEN
              IF (IERR.EQ.0) THEN
                CALL GSPMCI (4)
                CALL GPM (1,XVAL,YVAL)
              ELSE IF (IERR.EQ.1) THEN
                CALL GSPMCI (5)
                CALL GPM (1,XVAL,YVAL)
              ELSE
                CALL GSPMCI (6)
                CALL GPM (1,XVAL,YVAL)
              END IF
            ELSE
              IF (ABS(RLAT(I,J)).GT. 90..OR.ABS(RLON(I,J)).GT.180.) THEN
                RLAT(I,J)=1.E12
                RLON(I,J)=1.E12
                CALL GSPMCI (7)
                CALL GPM (1,XVAL,YVAL)
              END IF
            END IF
  101     CONTINUE
  102   CONTINUE
C
        CALL CPRECT (RLAT,MPTS,NPTS,NPTS,RWRK,LRWK,IWRK,LIWK)
        CALL CPCLDR (RLAT,RWRK,IWRK)
        CALL CPRECT (RLON,MPTS,NPTS,NPTS,RWRK,LRWK,IWRK,LIWK)
        CALL CPCLDR (RLON,RWRK,IWRK)
C
        CALL SFLUSH
        CALL GSPLCI (1)
        CALL GSTXCI (1)
        CALL GSMKSC (1.)
C
C Done.
C
        RETURN
C
      END


      SUBROUTINE MDRGDI (DINM)
C
C This is a user-replaceable routine that returns the name of the
C directory in which the RANGS/GSHHS data files have been placed.
C
        CHARACTER*(*) DINM
C
C Declare a common block in which the directory name will be passed.
C
        COMMON /EDRGDI/ USNM
          CHARACTER*60 USNM
        SAVE   /EDRGDI/
C
C Return the name of the directory.
C
        DINM=USNM
C
C Done.
C
        RETURN
C
      END


      FUNCTION IGETAI (IAMA,XCOP,YCOP,IAAI,IAGI,MNAI,NOAI,IGRP)
C
C This function retrieves from an area map the area identifier, relative
C to group IGRP, of the area containing the point (XCOP,YCOP).
C
        DIMENSION IAAI(MNAI),IAGI(MNAI)
C
        CALL ARGTAI (IAMA,XCOP,YCOP,IAAI,IAGI,MNAI,NOAI,1)
C
        DO 101 I=1,MNAI
          IF (IAGI(I).EQ.IGRP) THEN
            IGETAI=IAAI(I)
            RETURN
          END IF
  101   CONTINUE
C
        IGETAI=-1
        RETURN
C
      END
