      PARAMETER(RE = 6.8E8, PI=3.141592)
      PARAMETER(N=181, M=91, IMAX=4320, JMAX=2160, HPOLE=2810.)
      PARAMETER(IBR = ((N+2)*(M+2)*2/32+1)*3)
      PARAMETER(XSIZE=25., YSIZE=18.)
      PARAMETER(VXMIN=0.1, VXMAX=.9, VYMIN=0.05, VYMAX=0.65)
      PARAMETER(TXH=0.015)

      REAL       H(N,M)
      INTEGER*2 IH(IMAX)
      INTEGER   IW(IBR)
      CHARACTER CTRF*18, CSCALE*12, CTMP*24, CLAT*24, CLON*24
      EXTERNAL  IMOD

      CALL GLISET('MSGLEV', 2)

      WRITE(*, *) 'XLON, YLAT, ISCALE'
      READ (*, *)  XLON, YLAT, ISCALE
*      WRITE(*, *) 'SCALE (x 1,000,000)'
*      READ (*, *)  ISCALE

      CALL SCPARM(N, M, XLON, YLAT, ISCALE, 
     #            IX1, JY1, IDX, JDY, IMX, JMY)
*      WRITE(6,*) 'ORIGIN  ', IX1, JY1
*      WRITE(6,*) 'INTERVAL', IDX, JDY
*      WRITE(6,*) 'GRID NUM', IMX, JMY

      UXMIN = (IX1-1)/12.
      UXMAX = (IX1+IDX*(IMX-1)-1)/12.

      UYMAX = 90. - (JY1-1)/12.
      UYMIN = 90. - (JY1+JDY*(JMY-1)-1)/12.

      DX   = UXMAX - UXMIN
      DY   = UYMAX - UYMIN

      R = RE/(ISCALE*1.E6)/XSIZE

*------------------------ READING DATA ----------------------------------

      OPEN(10, FILE='etopo5.dat', ACCESS='DIRECT', RECL=8640)

      DO 10 J=1, JMY
        JREC = JY1 + (JMY-J)*JDY
        IF(1 .LE. JREC .AND. JREC .LE. JMAX) THEN
          READ(10, REC=JREC) IH
          DO 20 I=1, IMX
            IREC = IMOD(IX1+(I-1)*IDX-1, IMAX) + 1
            IF(1.LE.IREC .AND. IREC.LE.IMAX) THEN
              H(I,J) = IH(IREC)
            ELSE
              CALL MSGDMP('E', 'ETOPO', 'INVALID REC NO. (I)')
            ENDIF
   20     CONTINUE
        ELSE IF(JREC.EQ.JMAX+1) THEN
          DO 30 I=1, IMX
            H(I,J) = HPOLE
   30     CONTINUE
        ELSE
          CALL MSGDMP('E', 'ETOPO', 'INVALID REC NO. (J)')
        ENDIF
   10 CONTINUE

*------------------------ GRAPHIC OUTPUT ----------------------------------

      IF(ISCALE.GE.200) THEN
        ITR = 15
        PLY = 90.
        DH = 2000.
      ELSE
        ITR = 33
        PLY = YLAT
        DH  = 1000.
      ENDIF

      CALL SGPWSN
      READ (5,*) IWS
      CALL GROPN(IWS)

      CALL SGLSET('LFULL', .TRUE.)
      CALL SLFORM(XSIZE, YSIZE)
      CALL SLRAT (XSIZE, YSIZE)

      CALL GRFRM

      CALL SGSMPL(XLON, PLY, 0.)
      CALL SGSSIM(R, 0., 0.)
      CALL SGSTRN(ITR)
      CALL SGSVPT(VXMIN, VXMAX, VYMIN, VYMAX)
      CALL SGSWND(UXMIN, UXMAX, UYMIN, UYMAX)
      CALL SGSTXY(-180., 180., -90., 90.)
      CALL SGSTRF

      CALL SLPVPR(1)

*---------------------------- CONTOUR LINE -----------------------------

      CALL SGLSET('LCLIP', .TRUE.)
      CALL UDLSET('LABEL', .FALSE.)
      CALL UDLSET('LMSG',  .FALSE.)
      CALL UDLSET('LDASH', .FALSE.)

      CALL UMPGLB
      CALL UDGCLA(-10000., 10000., DH)
      CALL UDCNTZ(H, IMX, N, JMY, IW, IBR )

*----------------------------- MESSAGE --------------------------------

      CALL SGLSET('LCLIP', .FALSE.)

      CALL SGTRNL(ITR, CTRF)
      WRITE(CTMP, '(I3,A)') ISCALE, ',000,000'
      CALL CLADJ(CTMP)
      CSCALE = '1:'//CTMP
      CALL SGTXZV(VXMIN, VYMAX+    TXH, CSCALE, TXH, 0, -1, 1)
      CALL SGTXZV(VXMIN, VYMAX+2.5*TXH, CTRF,   TXH, 0, -1, 1)

      WRITE(CLON, '(A,F6.1)') 'CENTER LON.:', XLON
      WRITE(CLAT, '(A,F6.1)') '       LAT.:', YLAT
      CALL SGTXZV(VXMAX, VYMAX+    TXH, CLAT,   TXH, 0,  1, 1)
      CALL SGTXZV(VXMAX, VYMAX+2.5*TXH, CLON,   TXH, 0,  1, 1)

      WRITE(CTMP, '(A,F5.0,A)') 'CONTOUR INT.:', DH, 'm'
      CALL SGTXZV(VXMIN, VYMIN-    TXH, CTMP,   TXH, 0, -1, 1)

      CALL UMRGET('DGRIDMN', DGRID)
      WRITE(CTMP, '(A,F6.1)') 'GRID INT.:', DGRID
      CALL SGTXZV(VXMAX, VYMIN-    TXH, CTMP,   TXH, 0,  1, 1)


      CALL GRCLS
      END

*----------------------------------------------------------------------
      SUBROUTINE SCPARM(N, M, XLON, YLAT, ISCALE,
     #                  IX1, JY1, IDX, JDY, IMX, JMY)

      PARAMETER(PI=3.141592, IMAX=4320, JMAX=2160)

      RCT = COS(YLAT*PI/180.)
      IF(RCT.LE.0.05) THEN
        IFAC = 48
      ELSE IF(RCT.LE.0.1) THEN
        IFAC = 24
      ELSE IF(RCT.LE.0.25) THEN
        IFAC = 8
      ELSE IF(RCT.LE.0.34) THEN
        IFAC = 6
      ELSE IF(RCT.LE.0.5) THEN
        IFAC = 4
      ELSE IF(RCT.LE.0.67) THEN
        IFAC = 3
      ELSE
        IFAC = 2
      END IF

      IF(ISCALE.GE.200) THEN
        ISCALE = 200
      ELSE IF(ISCALE.GE.100) THEN
        ISCALE = 100
      ELSE IF(ISCALE.GE.50) THEN
        ISCALE = 50
      ELSE IF(ISCALE.GE.25) THEN
        ISCALE = 25
      ELSE IF(ISCALE.GE.10) THEN
        ISCALE = 10
      ELSE IF(ISCALE.GE.5) THEN
        ISCALE = 5
      ELSE
        ISCALE = 2
      END IF

      IDMAX = IMAX/(N-1)
      JDMAX = JMAX/(M-1)

      JDY  = MIN(JDMAX, MAX(1, JDMAX*ISCALE/90))
      IDX  = MIN(IDMAX, IFAC*JDY/2)

      ICT = NINT(XLON*12)      + 1
      JCT = NINT((90-YLAT)*12) + 1

      IX1 = ICT - (N-1)/2*IDX
      JY1 = JCT - (M-1)/2*JDY
      IF(JY1.LT.1) JY1 = 1

      IX2 = ICT + (N-1)/2*IDX
      JY2 = JCT + (M-1)/2*JDY
      IF(JY2.GT.JMAX) JY2 = JMAX

      IMX = MIN((IX2-IX1)/IDX + 1, IMAX/IDX+1)
      JMY = (JY2-JY1)/JDY + 1

      END
