PROGRAM CR2PLOT
C
C READ CR-2 MODE DETAIL VOYAGER MAG DATA AND PLOT
C
C ORIGINAL CODE - S. KRAMER  11/05/96
C
      PARAMETER (IARR=100000)
C
C PLOTTING PARAMETERS
C
      CHARACTER PLOTTIME*8,PLOTDATE*9,RUNTIME*12
      CHARACTER*10 VALUE,XFORM,YFORM,CTIME1,CTIME2
      CHARACTER*50 XSTRING,YSTRING,TITLE,SUBTITLE
      INTEGER*4 XLAB,YLAB,TYPE
      LOGICAL*1 DONE
      REAL*4 B(IARR),BX(IARR),BY(IARR),BZ(IARR),BDEL(IARR),BLAM(IARR),
     &       TIME(IARR),TDAY(IARR),THOUR(IARR)
C
C DATA INPUT PARAMETERS
C
      CHARACTER RECTYPE*4,TELFMT*4,FLTID*4,TIMEFMT*4,RUNMONTH*4,
     &          DSN*50,MODE*4,COORD*2
      INTEGER*2 DATE(6),DETAIL1(6),DETAIL2(6),TBEG(6),TEND(6),
     &          MODCNT(3),DATAID(2)
      INTEGER*4 RUNDAY,RUNYEAR
      LOGICAL*1 VOYSYS(32),MAGSYS(32),SYS1(32),SYS2(32),REC(10000)
      REAL*4 DETOUT(14432),DATA(14400),HDR(32),HDR1(100),
     &       PRIDAT(400,3),SECDAT(200,3),SCFDAT(400,3),AMBDAT(400,3),
     &       DATA1(IARR,3)
      REAL*8 DD,TD,TN,TP,DECDY,DECHR
C
      EQUIVALENCE (HDR(1),RECTYPE),    (HDR(2),TELFMT),
     &            (HDR(3),FLTID),      (HDR(4),DATE(1)),
     &            (HDR(7),DD),         (HDR(9),TD),
     &            (HDR(11),TIMEFMT),   (HDR(12),TIMEPD),
     &            (HDR(13),MODCNT(1)), (HDR(17),DATAID(1))
C
      EQUIVALENCE (HDR1(7),RUNYEAR),   (HDR1(9),RUNMONTH),
     &            (HDR1(10),RUNDAY),   (HDR1(69),SYS1(1)),
     &            (HDR1(77),SYS2(1)),  (DETOUT(1),HDR1(1))
C
      EQUIVALENCE (DETOUT(1),HDR(1)),     (DETOUT(33),DATA(1)),
     &            (PRIDAT(1,1),DATA(1)),  (SECDAT(1,1),DATA(1201)),
     &            (AMBDAT(1,1),DATA(1)),  (SCFDAT(1,1),DATA(1201))
C
      EQUIVALENCE (DETOUT(1),REC(1))
C
      DATA TITLE/' '/, SUBTITLE/' '/
C
C CR-2 PRIMARY DATA RESOLUTION
C
      DETAIL1(1) = 0 
      DETAIL1(2) = 0
      DETAIL1(3) = 0
      DETAIL1(4) = 0
      DETAIL1(5) = 0
      DETAIL1(6) = 120
C
C CR-2 SECONDARY DATA RESOLUTION
C
      DETAIL2(1) = 0 
      DETAIL2(2) = 0
      DETAIL2(3) = 0
      DETAIL2(4) = 0
      DETAIL2(5) = 0
      DETAIL2(6) = 240
C
      WRITE(6,*) '   ENTER DETAIL DSN'
      READ(5,'(A)') DSN
      OPEN(70,FILE=DSN,STATUS='OLD',FORM='FORMATTED',
     &     RECORDTYPE='VARIABLE',RECL=8191,READONLY)
      WRITE(6,'(1X,A)') DSN
C
      WRITE(6,*)
      WRITE(6,*) '   ENTER START TIME'
      WRITE(6,*) 'YY/DDD/HH:MM:SS.SSS'
      READ(5,'(I2,1X,I3,3(1X,I2),1X,I3)') TBEG
      WRITE(6,'(1X,I2,1X,I3,3(1X,I2),1X,I3)') TBEG
      WRITE(6,*)
      WRITE(6,*) '   ENTER STOP TIME'
      WRITE(6,*) 'YY/DDD/HH:MM:SS.SSS'
      READ(5,'(I2,1X,I3,3(1X,I2),1X,I3)') TEND
      WRITE(6,'(1X,I2,1X,I3,3(1X,I2),1X,I3)') TEND
C
      IDAYS = TEND(2) - TBEG(2)
      XL = DECHR(TBEG)
      XH = DECHR(TEND) + 24.0*FLOAT(IDAYS)
      WRITE(SUBTITLE,'(I4,'', DAY '',I3.3,A2)') 
     & TBEG(1)+1900,TBEG(2),'\E'
C
      FILL = 999.0
      NCNT = 0
      NPT1 = 0
      NPT2 = 0
   10 CONTINUE
       READ(70,'(Q,<LEN>A1)',END=100,ERR=10) LEN,(REC(I),I=1,LEN)
       NCNT = NCNT + 1
C
C EXTRACT HDR1 INFO
C
       IF ( DATAID(1).EQ.8 .AND. NCNT.EQ.1 ) THEN
C
        DO II = 1,32
         MAGSYS(II) = SYS2(II)
         VOYSYS(II) = SYS1(II)
C        WRITE(6,'(1X,''MAGSYS('',I2,'') = '',L1)') II,MAGSYS(II)
        END DO
C
        IF ( MAGSYS(4).OR.MAGSYS(5) ) THEN
         CONTINUE
        ELSE
         WRITE(6,*)
         WRITE(6,*) 'NOT A DETAIL DATA SET'
         STOP
        END IF
C
        IF ( MAGSYS(7) ) THEN
         COORD = 'HG'
         IF ( MAGSYS(17) ) COORD = 'S3'
         IF ( MAGSYS(18) ) COORD = 'L1'
         IF ( MAGSYS(19) ) COORD = 'U1'
         IF ( MAGSYS(20) ) COORD = 'N1'
         IF ( VOYSYS(8) ) COORD = 'HG'
         IF ( VOYSYS(9) ) COORD = 'L1'
         IF ( VOYSYS(10) ) COORD = 'S3'
         IF ( VOYSYS(11) ) COORD = 'U1'
         IF ( VOYSYS(12) ) COORD = 'N1'
        ELSE
         COORD = 'PL'
        END IF
C
        WRITE(6,*)
        WRITE(6,801) NCNT,RECTYPE,DATE,FLTID,COORD
C
        WRITE(RUNTIME,'(A4,''/'',I2.2,''/'',I4)')
     &   RUNMONTH,RUNDAY,RUNYEAR
C
C SELECT DESIRED DATA BLOCK
C
        WRITE(6,*)
        WRITE(6,*) 'ENTER PLOT DESIRED:'
        IF ( MAGSYS(4) ) THEN
         WRITE(6,*) '  1 - PRIMARY'
         WRITE(6,*) '  2 - SECONDARY'
        ELSE
         WRITE(6,*) '  3 - AMBIENT'
         WRITE(6,*) '  4 - S/C FIELD'
        END IF
        READ(5,*) IDAT
        IF ( IDAT.EQ.1 ) MODE = ' PRI'
        IF ( IDAT.EQ.2 ) MODE = ' SEC'
        IF ( IDAT.EQ.3 ) MODE = ' AMB'
        IF ( IDAT.EQ.4 ) MODE = ' SCF'
C
        READ(FLTID,'(3X,I1)') ID
        IF (FLTID.EQ.'FLT1') THEN
         TITLE = 'VOYAGER 1 '//COORD//MODE//'\E'
        ELSE IF (FLTID.EQ.'FLT2') THEN
         TITLE = 'VOYAGER 2 '//COORD//MODE//'\E'
        ELSE
         WRITE(6,*) '***INVALID FLIGHT ID***',FLTID,'***'
         STOP
        END IF
C
       END IF
C
C ASSUME TIME REQUEST << 1 YR
C
       IF ( DATE(1).LT.TBEG(1) ) GOTO 10
       IF ( DATE(2).LT.TBEG(2) ) GOTO 10
       IF ( DATE(2).EQ.TBEG(2) .AND.
     &      DECHR(DATE).LT.DECHR(TBEG) ) GOTO 10
       IF ( DATE(2).GE.TEND(2).AND.
     &      DECHR(DATE).GE.DECHR(TEND) ) GOTO 100       
C
       IF (DATAID(1).NE.1) GOTO 10
C
C PROCESS CR-2 LOW FIELD MAG DATA
C       
       DO I = 1,400
C
        ISEC = (I-1)/2 + 1
C
C GET DETAIL PRI/SEC DATA
C
        IF ( MAGSYS(4) ) THEN
         NPT1 = NPT1 + 1
         IF ( IDAT.EQ.1) THEN
          DATA1(NPT1,1) = PRIDAT(I,1)
          DATA1(NPT1,2) = PRIDAT(I,2)
          DATA1(NPT1,3) = PRIDAT(I,3)
         END IF
         IF ( IDAT.EQ.2 .AND. MOD(I-1,8).EQ.0 ) THEN
          NPT2 = NPT2 + 1
          TDAY(NPT2) = DECDY(DATE)
          THOUR(NPT2) = DECHR(DATE) + 24.0*(DATE(2)-TBEG(2))
          CALL INC_TIME(DATE,DETAIL2)
          DATA1(NPT2,1) = SECDAT(ISEC,1)
          DATA1(NPT2,2) = SECDAT(ISEC,2)
          DATA1(NPT2,3) = SECDAT(ISEC,3)
         END IF 
C
C GET DETAIL AMB/SCF DATA
C
        ELSE IF ( MAGSYS(5) ) THEN
C        WRITE(6,800) 
C    &    ID,DATE,I,(AMBDAT(I,J),J=1,3),I,(SCFDAT(I,J),J=1,3)
         NPT1 = NPT1 + 1
         IF (IDAT.EQ.3) THEN
          DATA1(NPT1,1) = AMBDAT(I,1)
          DATA1(NPT1,2) = AMBDAT(I,2)
          DATA1(NPT1,3) = AMBDAT(I,3)
         ELSE IF (IDAT.EQ.4) THEN
          DATA1(NPT1,1) = SCFDAT(I,1)
          DATA1(NPT1,2) = SCFDAT(I,2)
          DATA1(NPT1,3) = SCFDAT(I,3)
         END IF
        END IF
C
        IF ( IDAT.NE.2 ) THEN
         TDAY(NPT1) = DECDY(DATE)
         THOUR(NPT1) = DECHR(DATE) + 24.0*(DATE(2)-TBEG(2))
         CALL INC_TIME(DATE,DETAIL1)
        END IF
C
       END DO
C
      GOTO 10
  100 CONTINUE
C
      WRITE(6,*)
      WRITE(6,*) NCNT, ' RECORDS READ'
      WRITE(6,*) NPT1, ' PRI/AMB/SCF POINTS'
      WRITE(6,*) NPT2, ' SEC POINTS'
C
      IF ( NCNT.LE.0.OR.NCNT.GT.IARR ) GOTO 999
      IF ( IDAT.NE.2 .AND. (NPT1.LE.0.OR.NPT1.GT.IARR) ) GOTO 999
      IF ( IDAT.EQ.2 .AND. (NPT2.LE.0.OR.NPT2.GT.IARR) ) GOTO 999
C
      BMIN = 999.
      BXMIN = 999.
      BYMIN = 999.
      BZMIN = 999.
      BMAX = -999.
      BXMAX = -999.
      BYMAX = -999.
      BZMAX = -999.
C
C PRI/SEC/AMB/SCF PLOT DATA
C
      IF ( IDAT.NE.2 ) NMAX = NPT1
      IF ( IDAT.EQ.2 ) NMAX = NPT2
      DO I = 1,NMAX
C
C ASSIGN DECIMAL HOUR TO TIME ARRAY
C
       TIME(I) = THOUR(I)
C
C ASSIGN DECIMAL DAY TO TIME ARRAY
C
C       TIME(I) = TDAY(I)
C       WRITE(6,805) 
C    &   ID,TIME(I),I,(DATA1(I,J),J=1,3)
       BX(I) = DATA1(I,1)
       BY(I) = DATA1(I,2)
       BZ(I) = DATA1(I,3)
       IF ( BX(I).NE.FILL.AND.
     &      BY(I).NE.FILL.AND.
     &      BZ(I).NE.FILL) THEN
C
C GET MIN/MAX FIELD COMPONENTS
C
        IF ( BX(I).LT.BXMIN ) BXMIN = BX(I)
        IF ( BY(I).LT.BYMIN ) BYMIN = BY(I)
        IF ( BZ(I).LT.BZMIN ) BZMIN = BZ(I)
        IF ( BX(I).GT.BXMAX ) BXMAX = BX(I)
        IF ( BY(I).GT.BYMAX ) BYMAX = BY(I)
        IF ( BZ(I).GT.BZMAX ) BZMAX = BZ(I)
C
C COMPUTE FIELD MODULUS
C
        B2 = SQRT(BX(I)**2+BY(I)**2+BZ(I)**2)
        B(I) = B2
        IF ( B2.LT.BMIN ) BMIN = B2
        IF ( B2.GT.BMAX ) BMAX = B2
C
C COMPUTE LATITUDINAL ANGLE
C
        IF (B2.NE.0.0) THEN
         BDEL(I) = ASIN(BZ(I)/B2) * 180.0/3.14159
        ELSE
         BDEL(I) = FILL
        END IF
C
C COMPUTE LONGITUDINAL ANGLE
C
        IF (BX(I).NE.0.0.OR.BY(I).NE.0.0) THEN
         BLAM(I) = 180.0 - ATAN2(BY(I),-BX(I))  * 180.0/3.14159
        END IF
       ELSE
        B(I) = FILL
        BDEL(I) = FILL
        BLAM(I) = FILL
       END IF
C
C END PRI/SEC/AMB/SCF PLOT DATA LOOP
C
      END DO
C
      WRITE(6,*)
      WRITE(6,'(1X,''NUMBER OF RECORDS READ: '',I5)') NCNT
      WRITE(6,'(1X,''NUMBER OF POINTS EXTRACTED: '',I5)') NMAX
      CLOSE(70)
C
  800 FORMAT(1X,I1,1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3,
     &       2(1X,I6,3(1X,F7.3)))
  801 FORMAT(1X,I5,1X,A4,1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3,1X,A4,1X,A2)
  802 FORMAT(1X,I2,1X,I3.3,3(1X,I2.2),1X,I3.3,2(1X,I4,3(1X,F7.3)))
  805 FORMAT(1X,I1,1X,F12.8,1X,I6,3(1X,F7.3))
C
C MAG PLOT ROUTINE
C
      NPTS = NMAX
      WRITE(6,*) 
      WRITE(6,*) NMAX,' POINTS PER PLOT'
C
      CALL GETTIME(PLOTTIME)
      CALL GETDATE(PLOTDATE)
C
      WRITE(6,*)
      WRITE(6,*) 'ENTER PLOT DISCONTINUITY INTERVAL (HOURS)'
      READ(5,*) DTIME
C
      IF (TITLE(1:1).EQ.' ') THEN
       WRITE(6,*)
       WRITE(6,*) 'ENTER TITLE'
       READ(5,'(Q,A)') LEN,TITLE
       TITLE(LEN+1:) = '\E'
      END IF
      IF (SUBTITLE(1:1).EQ.' ') THEN
       WRITE(6,*)
       WRITE(6,*) 'ENTER SUBTITLE'
       READ(5,'(Q,A)') LEN,SUBTITLE
       SUBTITLE(LEN+1:) = '\E'
      END IF
C
      CALL MGOINIT
      CALL MGOSETUP(-5)
      CALL MGOERASE
C
      XX1 = 400.
      XX2 = 2225.
C
C PLOT SUBTITLE
C
      YY1 = 3000.
      YY2 = 3100.
      CALL MGOSETEXPAND(1.0)
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(0.,0.,100.,100.)
      CALL MGORELOCATE(50.,50.)
      CALL MGOPUTLABEL(50,SUBTITLE,5)
C
C FIELD MAGNITUDE
C
      YY1 = 2730.
      YY2 = 3000.
      XT1 = 100.
      XT2 = 200.
      YT1 = 100.
      YT2 = 200.
      TT1 = 3100.
      TT2 = 3200.
      XNUMH = .8
      YNUMH = .8
      XTXTH = 1.
      YTXTH = 1.
      TITLH = 1.3
      XLAB = 4
      YLAB = 3
      XTLEN = 50.
      YTLEN = 50.
      TYPE = 0
C
      XSTRING = 'HOUR\E'
C     XSTRING = 'DAY\E'
      XLOW = XL
      XHIGH = XH
      WRITE(6,*)
      WRITE(6,*) 'ENTER BIG TICK SPACING   (DECIMAL HOURS)'
      READ(5,*) BSPANX
      WRITE(6,*)
      WRITE(6,*) 'ENTER SMALL TICK SPACING (DECIMAL MINUTES)'
      READ(5,*) SSPANX
      SSPANX = SSPANX/60.0      
      XFORM = 'F5.1'
C     WRITE(6,*) XLOW,XHIGH,BSPANX,SSPANX
C     WRITE(6,*) TIME(1),TIME(NPTS),NPTS
      YSTRING = 'B (nT)\E'
C     YLOW = 0.
C     YHIGH = .6
C     BSPANY = .2
C     SSPANY = .05
      WRITE(6,*)
      WRITE(6,'(1X,''BMIN:  '',F7.3,3X,''BMAX:  '',F7.3)') 
     & BMIN,BMAX
      WRITE(6,*)
      WRITE(6,*) 'ENTER FIELD MAGNITUDE LIMITS'
      WRITE(6,*) '(YLOW,YHIGH,BSPANY,SSPANY)'
      READ(5,*) YLOW,YHIGH,BSPANY,SSPANY
      YFORM = 'F6.1'
C
C FLAG OUT OF BOUNDS POINTS AND REMOVE
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
      CALL MGOSETEXPAND(1.0)
      DO I = 1,NPTS
       IF (B(I).GT.YHIGH.AND.B(I).NE.FILL) THEN
        FOLD = MOD(B(I),YHIGH)
        CALL MGORELOCATE(TIME(I),FOLD)
        CALL MGOPOINT(10,1)
        B(I) = FILL
       END IF
      END DO
C
C PLOT GRID AND LABELS
C
      CALL GRID(XX1,XX2,YY1,YY2,XT1,XT2,YT1,YT2,TT1,TT2,
     C          XNUMH,YNUMH,XTXTH,YTXTH,TITLH,XLAB,YLAB,
     C          XTLEN,YTLEN,XSTRING,YSTRING,TITLE,TYPE,
     C          XLOW,XHIGH,YLOW,YHIGH,SSPANX,BSPANX,
     C          SSPANY,BSPANY,XFORM,YFORM)
C
C PLOT DATA
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
C     CALL MGOPLOTID(F1,F2)
      CALL MGOSETEXPAND(0.5)
      CALL MGOGRELOCATE(XX2,YY2)
      CALL MGOPUTLABEL(20,PLOTDATE//'   '//PLOTTIME,7)
      CALL MGOGRELOCATE(XX1,YY2)
      CALL MGOPUTLABEL(31,'EDR PROCESS DATE:  '//RUNTIME,9)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
C     CALL MGOCONNECT(TIME,B,NPTS)
      CALL MGOSETEXPAND(0.01)
C
      CALL MCONNECT(FILL,TIME,B,NPTS,DTIME)
C
C LATITUDINAL FIELD ANGLE
C
      YY1 = 2460.
      YY2 = 2730.
      XLAB = 4
      TITLH = 0.
      YSTRING = '\Gd\DB\E'
      YLOW = -90.
      YHIGH = 90.
      SSPANY = 15.
      BSPANY = 45.
C
C PLOT GRID AND LABELS
C
      CALL GRID(XX1,XX2,YY1,YY2,XT1,XT2,YT1,YT2,TT1,TT2,
     C          XNUMH,YNUMH,XTXTH,YTXTH,TITLH,XLAB,YLAB,
     C          XTLEN,YTLEN,XSTRING,YSTRING,TITLE,TYPE,
     C          XLOW,XHIGH,YLOW,YHIGH,SSPANX,BSPANX,
     C          SSPANY,BSPANY,XFORM,YFORM)
C
C PLOT DATA
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
C     CALL MGOCONNECT(TIME,BDEL,NPTS)
      CALL MGOSETEXPAND(0.01)
      CALL MCONNECT(FILL,TIME,BDEL,NPTS,DTIME)
      CALL MGOSETEXPAND(1.0)
C
C LONGITUDINAL FIELD ANGLE
C
      YY1 = 1920.
      YY2 = 2460.
      XLAB = 4
      TITLH = 0.
      YSTRING = '\Gl\DB\E'
      YLOW = 0.
      YHIGH = 360.
      SSPANY = 15.
      BSPANY = 90.
C
C PLOT GRID AND LABELS
C
      CALL GRID(XX1,XX2,YY1,YY2,XT1,XT2,YT1,YT2,TT1,TT2,
     C          XNUMH,YNUMH,XTXTH,YTXTH,TITLH,XLAB,YLAB,
     C          XTLEN,YTLEN,XSTRING,YSTRING,TITLE,TYPE,
     C          XLOW,XHIGH,YLOW,YHIGH,SSPANX,BSPANX,
     C          SSPANY,BSPANY,XFORM,YFORM)
C
C PLOT DATA
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
C     CALL MGOCONNECT(TIME,BLAM,NPTS)
      CALL MGOSETEXPAND(0.01)
      CALL MCONNECT(FILL,TIME,BLAM,NPTS,DTIME)
      CALL MGOSETEXPAND(1.0)
C
C X FIELD COMPONENT
C
      YY1 = 1380.
      YY2 = 1920.
      XLAB = 4
      TITLH = 0.
      IF ( COORD.EQ.'PL' ) THEN
       YSTRING = 'B\DX (nT)\E'
      ELSE IF ( COORD.EQ.'HG' ) THEN
       YSTRING = 'B\DR (nT)\E'
      ELSE
       YSTRING = 'B\DR (nT)\E'
      END IF
C     YLOW = -.6
C     YHIGH = .6
C     SSPANY = .05
C     BSPANY = .2
      WRITE(6,*)
      WRITE(6,'(1X,''BXMIN:  '',F8.3,3X,''BXMAX:  '',F8.3)') 
     & BXMIN,BXMAX
      WRITE(6,*)
      WRITE(6,*) 'ENTER X FIELD COMPONENT LIMITS'
      WRITE(6,*) '(YLOW,YHIGH,BSPANY,SSPANY)'
      READ(5,*) YLOW,YHIGH,BSPANY,SSPANY
C
C FLAG OUT OF BOUNDS POINTS AND REMOVE
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
      CALL MGOSETEXPAND(1.0)
      DO I = 1,NPTS
       IF ((BX(I).LT.YLOW.OR.
     &      BX(I).GT.YHIGH).AND.
     &      BX(I).NE.FILL) THEN
        IF (BX(I).LT.YLOW) FOLD = MOD(BX(I),YLOW)
        IF (BX(I).GT.YHIGH) FOLD = MOD(BX(I),YHIGH)
        CALL MGORELOCATE(TIME(I),FOLD)
        CALL MGOPOINT(10,1)
        BX(I) = FILL
       END IF
      END DO
C
C PLOT GRID AND LABELS
C
      CALL GRID(XX1,XX2,YY1,YY2,XT1,XT2,YT1,YT2,TT1,TT2,
     C          XNUMH,YNUMH,XTXTH,YTXTH,TITLH,XLAB,YLAB,
     C          XTLEN,YTLEN,XSTRING,YSTRING,TITLE,TYPE,
     C          XLOW,XHIGH,YLOW,YHIGH,SSPANX,BSPANX,
     C          SSPANY,BSPANY,XFORM,YFORM)
C
C PLOT DATA
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
C     CALL MGOCONNECT(TIME,BX,NPTS)
      CALL MGOSETEXPAND(0.01)
      CALL MCONNECT(FILL,TIME,BX,NPTS,DTIME)
C
C Y FIELD COMPONENT
C
      YY1 = 840.
      YY2 = 1380.
      XLAB = 4
      TITLH = 0.
      IF ( COORD.EQ.'PL' ) THEN
       YSTRING = 'B\DY (nT)\E'
      ELSE IF ( COORD.EQ.'HG' ) THEN
       YSTRING = 'B\DT (nT)\E'
      ELSE
       YSTRING = 'B\D\GQ (nT)\E'
      END IF
C     YLOW = -.6
C     YHIGH = .6
C     SSPANY = .05
C     BSPANY = .2
      WRITE(6,*)
      WRITE(6,'(1X,''BYMIN:  '',F8.3,3X,''BYMAX:  '',F7.3)') 
     & BYMIN,BYMAX
      WRITE(6,*)
      WRITE(6,*) 'ENTER Y FIELD COMPONENT LIMITS'
      WRITE(6,*) '(YLOW,YHIGH,BSPANY,SSPANY)'
      READ(5,*) YLOW,YHIGH,BSPANY,SSPANY
C
C FLAG OUT OF BOUNDS POINTS AND REMOVE
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
      CALL MGOSETEXPAND(1.0)
      DO I = 1,NPTS
       IF ((BY(I).LT.YLOW.OR.
     &      BY(I).GT.YHIGH).AND.
     &      BY(I).NE.FILL) THEN
        IF (BY(I).LT.YLOW) FOLD = MOD(BY(I),YLOW)
        IF (BY(I).GT.YHIGH) FOLD = MOD(BY(I),YHIGH)
        CALL MGORELOCATE(TIME(I),FOLD)
        CALL MGOPOINT(10,1)
        BY(I) = FILL
       END IF
      END DO
C
C PLOT GRID AND LABELS
C
      CALL GRID(XX1,XX2,YY1,YY2,XT1,XT2,YT1,YT2,TT1,TT2,
     C          XNUMH,YNUMH,XTXTH,YTXTH,TITLH,XLAB,YLAB,
     C          XTLEN,YTLEN,XSTRING,YSTRING,TITLE,TYPE,
     C          XLOW,XHIGH,YLOW,YHIGH,SSPANX,BSPANX,
     C          SSPANY,BSPANY,XFORM,YFORM)
C
C PLOT DATA
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
C     CALL MGOCONNECT(TIME,BY,NPTS)
      CALL MGOSETEXPAND(0.01)
      CALL MCONNECT(FILL,TIME,BY,NPTS,DTIME)
C
C Z FIELD COMPONENT
C
      YY1 = 300.
      YY2 = 840.
      XLAB = 0
      TITLH = 0.
      IF ( COORD.EQ.'PL' ) THEN
       YSTRING = 'B\DZ (nT)\E'
      ELSE IF ( COORD.EQ.'HG' ) THEN
       YSTRING = 'B\DN (nT)\E'
      ELSE
       YSTRING = 'B\D\GF (nT)\E'
      END IF
C     YLOW = -.6
C     YHIGH = .6
C     SSPANY = .05
C     BSPANY = .2
      WRITE(6,*)
      WRITE(6,'(1X,''BZMIN:  '',F8.3,3X,''BZMAX:  '',F7.3)') 
     & BZMIN,BZMAX
      WRITE(6,*)
      WRITE(6,*) 'ENTER Z FIELD COMPONENT LIMITS'
      WRITE(6,*) '(YLOW,YHIGH,BSPANY,SSPANY)'
      READ(5,*) YLOW,YHIGH,BSPANY,SSPANY
C
C FLAG OUT OF BOUNDS POINTS AND REMOVE
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
      CALL MGOSETEXPAND(1.0)
      DO I = 1,NPTS
       IF ((BZ(I).LT.YLOW.OR.
     &      BZ(I).GT.YHIGH).AND.
     &      BZ(I).NE.FILL) THEN
        IF (BZ(I).LT.YLOW) FOLD = MOD(BZ(I),YLOW)
        IF (BZ(I).GT.YHIGH) FOLD = MOD(BZ(I),YHIGH)
        CALL MGORELOCATE(TIME(I),FOLD)
        CALL MGOPOINT(10,1)
        BZ(I) = FILL
       END IF
      END DO
C
C PLOT GRID AND LABELS
C
      CALL GRID(XX1,XX2,YY1,YY2,XT1,XT2,YT1,YT2,TT1,TT2,
     C          XNUMH,YNUMH,XTXTH,YTXTH,TITLH,XLAB,YLAB,
     C          XTLEN,YTLEN,XSTRING,YSTRING,TITLE,TYPE,
     C          XLOW,XHIGH,YLOW,YHIGH,SSPANX,BSPANX,
     C          SSPANY,BSPANY,XFORM,YFORM)
C
C PLOT DATA
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETLIM(XL,YLOW,XH,YHIGH)
C     CALL MGOCONNECT(TIME,BZ,NPTS)
      CALL MGOSETEXPAND(0.01)
      CALL MCONNECT(FILL,TIME,BZ,NPTS,DTIME)
C
  888 CONTINUE
C
      CALL MGOPRNTPLOT(NVEC)
      WRITE(6,*)
      WRITE(6,'(1X,''VECTORS PLOTTED: '',I6)') NVEC
      WRITE(6,*)
C
  999 CONTINUE
C
      STOP
      END
      REAL*8 FUNCTION DECDY(DATE)
C
C CONVERT HIGH RESOLUTION CALENDAR DATA TO DECIMAL DAY
C
      INTEGER*2 DATE(6)
      REAL*8 DAY,HOUR,MIN,SEC,MSEC
      MSEC = DFLOAT(DATE(6))
      SEC = DFLOAT(DATE(5)) + MSEC/1000.0D0
      MIN = DFLOAT(DATE(4)) + SEC/60.0D0
      HOUR = DFLOAT(DATE(3)) + MIN/60.0D0
      DAY = DFLOAT(DATE(2)) + HOUR/24.0D0
      DECDY = DAY
      RETURN
      END
      REAL*8 FUNCTION DECHR(DATE)
C
C CONVERT HIGH RESOLUTION CALENDAR DATA TO DECIMAL HOUR
C
      INTEGER*2 DATE(6)
      REAL*8 DAY,HOUR,MIN,SEC,MSEC
      MSEC = DFLOAT(DATE(6))
      SEC = DFLOAT(DATE(5)) + MSEC/1000.0D0
      MIN = DFLOAT(DATE(4)) + SEC/60.0D0
      HOUR = DFLOAT(DATE(3)) + MIN/60.0D0
      DECHR = HOUR
      RETURN
      END
      SUBROUTINE INC_TIME(TIME,DELTA)
C
C INCREMENT CALENDAR TIME (YY,DDD,HH,MM,SS,FFF) BY DELTA(6)
C
      INTEGER*2 TIME(6),DELTA(6)
C
      IYR = TIME(1)
      LEAP = 365
      IF (MOD(IYR,4).EQ.0) LEAP = 366
C
      TIME(6) = TIME(6) + DELTA(6)
      IF (TIME(6).GT.999) THEN
       TIME(5) = TIME(5) + 1
       TIME(6) = TIME(6) - 1000
      END IF
C
      TIME(5) = TIME(5) + DELTA(5)
      IF (TIME(5).GT.59) THEN
       TIME(4) = TIME (4) + 1
       TIME(5) = TIME (5) - 60
      END IF
C
      TIME(4) = TIME(4) + DELTA(4)
      IF (TIME(4).GT.59) THEN
       TIME(3) = TIME(3) + 1
       TIME(4) = TIME(4) - 60
      END IF
C
      TIME(3) = TIME(3) + DELTA(3)
      IF (TIME(3).GT.23) THEN
       TIME(2) = TIME(2) + 1
       TIME(3) = TIME(3) - 24
      END IF
C
      TIME(2) = TIME(2) + DELTA(2)
      IF (TIME(2).GT.LEAP) THEN
       TIME(1) = TIME(1) + 1
       TIME(2) = TIME(2) - LEAP
       IF (TIME(1).GT.99) TIME(1) = TIME(1) - 100
      END IF
C
      TIME(1) = TIME(1) + DELTA(1)
C
      RETURN
      END
      REAL*8 FUNCTION DECYR(DATE)
C
C CONVERT HIGH RESOLUTION CALENDAR TIME INTO DECIMAL TIME
C
      INTEGER*2 DATE(6)
      REAL*8 YEAR,DAY,HOUR,MIN,SEC,MSEC,LEAP
      LEAP = 365.0D0
      IF (MOD(DATE(1),4).EQ.0) LEAP = 366.0D0
      MSEC = DFLOAT(DATE(6))
      SEC = DFLOAT(DATE(5)) + MSEC/1000.0D0
      MIN = DFLOAT(DATE(4)) + SEC/60.0D0
      HOUR = DFLOAT(DATE(3)) + MIN/60.0D0
      DAY = DFLOAT(DATE(2)-1) + HOUR/24.0D0
      YEAR = DFLOAT(DATE(1)) + DAY/LEAP
      DECYR = YEAR
      RETURN
      END
      REAL*4 FUNCTION DVAL(VAL)
C
C RETURN REAL*4 DECIMAL DAY FROM CALENDAR INPUT (YYDDDHH)
C
      CHARACTER*10 VALUE
      IVAL = INT(VAL)
      WRITE(VALUE,'(I7)') IVAL
      READ(VALUE,'(I2,I3.3,I2.2)') IY,ID,IH
      DVAL = FLOAT(ID) + FLOAT(IH)/24.0
      RETURN
      END
C
C MONGO PLOTTING ROUTINE WITH USER DEFINED SPECIFICS
C DEVELOPED BY SAUNDERS B. KRAMER JR., CODE 692
C
C XX1     PLOT LOWER X AXIS LIMIT IN DEVICE COORDINATES - R*4
C XX2     PLOT UPPER X AXIS LIMIT IN DEVICE COORDINATES - R*4
C YY1     PLOT LOWER Y AXIS LIMIT IN DEVICE COORDINATES - R*4
C YY2     PLOT UPPER Y AXIS LIMIT IN DEVICE COORDINATES - R*4
C XT1     Y LABEL LOWER X AXIS LIMIT IN DEVICE COORDINATES - R*4
C XT2     Y LABEL UPPER X AXIS LIMIT IN DEVICE COORDINATES - R*4
C YT1     X LABEL LOWER Y AXIS LIMIT IN DEVICE COORDINATES - R*4
C YT2     X LABEL UPPER Y AXIS LIMIT IN DEVICE COORDINATES - R*4
C TT1     TITLE LOWER Y AXIS LIMIT IN DEVICE COORDINATES - R*4
C TT2     TITLE UPPER Y AXIS LIMIT IN DEVICE COORDINATES - R*4
C XNUMH   SIZE OF X SCALE CHARACTERS - R*4
C YNUMH   SIZE OF Y SCALE CHARACTERS - R*4
C XTXTH   SIZE OF X LABEL CHARACTERS - R*4
C YTXTH   SIZE OF Y LABEL CHARACTERS - R*4
C TITLH   SIZE OF TITLE CHARACTERS - R*4
C XLAB    CONTROLS X AXIS SCALING AND LABELLING - I*4
C   0 PRINT ALL SCALES AND X LABEL
C   1 SKIP FIRST SCALE
C   2 SKIP LAST SCALE
C   3 SKIP FIRST AND LAST SCALE
C   4 NO SCALES 
C YLAB    PARALLELS XLAB ABOVE FOR Y AXIS - I*4
C XTLEN   SIZE OF X AXIS LARGE TICKS - R*4
C YTLEN   SIZE OF Y AXIS LARGE TICKS - R*4
C XSTRING X LABEL STRING - C*50
C YSTRING Y LABEL STRING - C*50
C TITLE TITLE STRING - C*50
C TYPE INDICATES INDEPENDENT AXIS FORMAT - I*4
C  -6 TIME FORMAT (YEAR//DAY) YEAR INTERVALS
C  -5 TIME FORMAT (YEAR)       *
C  -4 TIME FORMAT (HOUR)      DAY INTERVALS
C  -3 TIME FORMAT (DAY)        *  
C  -2 TIME FORMAT (DAY//HOUR)  *
C  -1 TIME FORMAT (YEAR//DAY)  *
C   0 STANDARD
C   1 LOGARITHMIC INDEPENDENT AXIS, STANDARD Y AXIS
C   2 LOGARITHMIC DEPENDENT AXIS, STANDARD X AXIS
C   3 LOGARITHMIC INDEPENDENT AND DEPENDENT AXIS
C XLOW    PLOT LOWER X AXIS LIMIT IN USER COORDINATES - R*4
C XHIGH   PLOT UPPER X AXIS LIMIT IN USER COORDINATES - R*4
C YLOW    PLOT LOWER Y AXIS LIMIT IN USER COORDINATES - R*4
C YHIGH   PLOT LOWER Y AXIS LIMIT IN USER COORDINATES - R*4
C SSPANX  SPACING OF INDEPENDENT AXIS SMALL TICKS - R*4
C BSPANX  SPACING OF INDEPENDENT AXIS LARGE TICKS - R*4
C SSPANY  SPACING OF DEPENDENT AXIS SMALL TICKS - R*4
C BSPANY  SPACING OF DEPENDENT AXIS LARGE TICKS - R*4
C XFORM   FORTRAN FORMAT OF INDEPENDENT AXIS SCALING - C*10
C YFORM   FORTRAN FORMAT OF DEPENDENT AXIS SCALING - C*10
C
      SUBROUTINE GRID(XX1,XX2,YY1,YY2,XT1,XT2,YT1,YT2,TT1,TT2,
     C                XNUMH,YNUMH,XTXTH,YTXTH,TITLH,XLAB,YLAB,
     C                XTLEN,YTLEN,XSTRING,YSTRING,TITLE,TYPE,
     C                XLOW,XHIGH,YLOW,YHIGH,SSPANX,BSPANX,
     C                SSPANY,BSPANY,XFORM,YFORM)
      CHARACTER*1 TCHAR
      CHARACTER*10 VALUE,XFORM,YFORM
      CHARACTER*50 XSTRING,YSTRING,TITLE
      INTEGER*4 XNUM,YNUM,XLAB,YLAB,TYPE
      LOGICAL TIME
      REAL*8 TLOW,THIGH,DEC,CONVERT
C
C DETERMINE IF INDEPENDENT AXIS HAS TIME SCALE
C
      TIME = .FALSE.
      IF (TYPE.LT.0) TIME = .TRUE.
      IF (.NOT.TIME) GOTO 30
C
C CONVERT CALENDAR TIME INTO DECIMAL YEAR
C
       TLOW = CONVERT(XLOW)
       THIGH = CONVERT(XHIGH)
   30 CONTINUE
C
C DRAW BOX
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOLINE(XX1,YY1,XX1,YY2)
      CALL MGOLINE(XX1,YY2,XX2,YY2)
      CALL MGOLINE(XX2,YY2,XX2,YY1)
      CALL MGOLINE(XX2,YY1,XX1,YY1)
C
C DRAW TIME (INDEPENDENT) AXIS TICK MARKS
C
      IF (.NOT.TIME) GOTO 40
       CALL MGOSETEXPAND(XNUMH)
       CALL BIGTICK(XX1,XX2,YY1,YY2,TLOW,THIGH,BSPANX,
     C              XLAB,XTLEN,TYPE)
       CALL SMALLTICK(XX1,XX2,YY1,YY2,TLOW,THIGH,SSPANX,
     C              XTLEN,TYPE)
      GOTO 45
   40 CONTINUE       
C
C DRAW INDEPENDENT AXIS LARGE TICK MARKS AND LABEL
C
      IF (TYPE.EQ.1.OR.TYPE.EQ.3) THEN
       ISCALE = 1
      ELSE 
       ISCALE = 0
      END IF
      CALL MGOSETEXPAND(XNUMH)
      CALL LRGTICK(YY1,YY2,XX1,XX2,XLOW,XHIGH,BSPANX,
     C             XFORM,XTLEN,XLAB,ISCALE)
C
C DRAW INDEPENDENT AXIS SMALL TICK MARKS
C
      IF (TYPE.EQ.1.OR.TYPE.EQ.3) THEN
       CALL LOGTICK(YY1,YY2,XX1,XX2,XLOW,XHIGH,SSPANX,XTLEN)
      ELSE
       CALL SMLTICK(YY1,YY2,XX1,XX2,XLOW,XHIGH,SSPANX,XTLEN)
      END IF
C
   45 CONTINUE
C 
C DRAW DEPENDENT AXIS LARGE TICK MARKS AND LABEL 
C
      IF (TYPE.EQ.2.OR.TYPE.EQ.3) THEN
       ISCALE = 1
      ELSE
       ISCALE = 0
      END IF
      CALL MGOSETEXPAND(YNUMH)
      XX1 = -XX1
      CALL LRGTICK(XX1,XX2,YY1,YY2,YLOW,YHIGH,
     C             BSPANY,YFORM,YTLEN,YLAB,ISCALE)
C
C DRAW DEPENDENT AXIS SMALL TICK MARKS 
C
      XX1 = -XX1
      IF (TYPE.EQ.2.OR.TYPE.EQ.3) THEN
       CALL LOGTICK(XX1,XX2,YY1,YY2,YLOW,YHIGH,SSPANY,YTLEN)
      ELSE
       CALL SMLTICK(XX1,XX2,YY1,YY2,YLOW,YHIGH,SSPANY,YTLEN)
      END IF
C
C LABEL TITLE
C
      CALL STRINGSET(XX1,TT1,XX2,TT2,TITLE,50,1,TITLH)
C
C LABEL INDEPENDENT AXIS
C
      CALL STRINGSET(XX1,YT1,XX2,YT2,XSTRING,50,1,XTXTH)
C
C LABEL DEPENDENT AXIS
C
      CALL MGOSETANGLE(90.)
      CALL STRINGSET(XT1,YY1,XT2,YY2,YSTRING,50,1,YTXTH)
      CALL MGOSETANGLE(0.)
C
C RESET LOCATION,SIZE AND WEIGHT
C
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      CALL MGOSETEXPAND(1.0)
      CALL MGOSETLWEIGHT(1)
C
      RETURN
      END
      SUBROUTINE STRINGSET(XX1,YY1,XX2,YY2,STRING,LEN,WGHT,SIZE)
C
C PLACE LABEL CENTERED IN A BOX DEFINED BY DEVICE COORDINATES
C
      CHARACTER*(*) STRING
      INTEGER*4 WGHT
      CALL MGOSETLOC(XX1,YY1,XX2,YY2)
      XX = (XX1 + XX2) / 2.0
      YY = (YY1 + YY2) / 2.0
      CALL MGOGRELOCATE(XX,YY)
      CALL MGOSETEXPAND(SIZE)
      CALL MGOSETLWEIGHT(WGHT)
      IF (SIZE.NE.0.) CALL MGOPUTLABEL(LEN,STRING,5)
      CALL MGOSETEXPAND(1.0)
      CALL MGOSETLWEIGHT(1)
      RETURN
      END
      SUBROUTINE BIGTICK(XX1,XX2,YY1,YY2,XLOW,XHIGH,SPAN,
     C                   XLAB,XTLEN,TYPE)
C
C PRODUCE LARGE TICK MARKS AND LABELS ALONG INDEPENDENT AXIS (TIME)
C AT USER DEFINED INTERVALS IN DAYS.  RANGE OF AXIS VALUES AND USER 
C SPECIFIED SPACING DETERMINE NUMBER OF TICK MARKS.  DEVICE 
C COORDINATES USED FOR LOCATING MARKS.
C
      CHARACTER*10 XCHAR
      INTEGER*4 XLAB,TYPE
      REAL*8 XLOW,XHIGH,XVAL,VAL,YR1,YR2,FRAC1,FRAC2
C
      Y0 = YY1 - XTLEN/2.
      Y1 = YY1 + XTLEN
      Y2 = YY2 - XTLEN
      IF (TYPE.GT.-5) THEN
       XDIFF = REAL((XHIGH-XLOW)*365.25D0)
      ELSE
       XDIFF = REAL(XHIGH-XLOW)
      END IF
      YDIFF = YY2 - YY1
      ISTOP = NINT(XDIFF/SPAN) + 1
      XVAL = XLOW
      DO 50 I=1,ISTOP
       XLIN = XX1 + (XX2-XX1) * REAL((XVAL-XLOW)/(XHIGH-XLOW))
       ERR = XVAL/XHIGH
       IF (ERR.GT.1.0) GOTO 50
C       ERR = REAL(XVAL-XHIGH)/REAL(XHIGH-XLOW)
C       IF (ERR.GT.0.0001) GOTO 50
C       IF (XVAL.GT.XHIGH) GOTO 50
       CALL MGOLINE(XLIN,YY1,XLIN,Y1)
       CALL MGOLINE(XLIN,YY2,XLIN,Y2)
       IF ( XLAB.EQ.4 ) GOTO 30
       CALL CALENDAR(XVAL,IY,ID,IH)
C BUILD SCALE STRING
       IF (TYPE.EQ.-1) WRITE(XCHAR,'(I2,I3.3)') IY,ID
       IF (TYPE.EQ.-2) WRITE(XCHAR,'(I3.3,I2.2)') ID,IH
       IF (TYPE.EQ.-3) WRITE(XCHAR,'(I3.3)') ID
       IF (TYPE.EQ.-4) WRITE(XCHAR,'(I2.2)') IH
       IF (TYPE.EQ.-5) WRITE(XCHAR,'(I2.2)') IY
       IF (TYPE.EQ.-6) WRITE(XCHAR,'(I2.2,I3.3)') IY,ID
       INUM = 0
       DO ILET=1,10
        LVAL = ICHAR(XCHAR(ILET:ILET))
        IF ( LVAL.NE.32 ) INUM = INUM + 1
       END DO
       IF ( XLAB.EQ.1.AND.I.EQ.1 ) GOTO 30
       IF ( XLAB.EQ.2.AND.I.EQ.ISTOP ) GOTO 30
       IF ( XLAB.EQ.3.AND.(I.EQ.1.OR.I.EQ.ISTOP) ) GOTO 30
       CALL MGOGRELOCATE(XLIN,Y0)
C       TYPE*,'X LABEL *'//XCHAR//'*'
       CALL MGOPUTLABEL(INUM,XCHAR,2)
   30  IF (I.EQ.ISTOP) GOTO 50
       IF (TYPE.GT.-5) THEN
C BALANCE CALENDAR ARITHMETIC ACROSS CHANGE OF YEAR
        YR1 = DBLE(INT(XVAL))
        VAL = XVAL + DBLE(SPAN) / DBLE(DAYS(XVAL))
        IF (INT(VAL).GT.INT(XVAL)) THEN
         YR2 = YR1 + 1.D0
         FRAC1 = (YR2 - XVAL) * DBLE(DAYS(XVAL))
         FRAC2 = DBLE(SPAN) - FRAC1 
         VAL = YR2 + FRAC2 / DBLE(DAYS(YR2))
        ENDIF
        XVAL = VAL
       ELSE
        XVAL = XVAL + SPAN
       END IF
   50 CONTINUE
C
      RETURN
      END
      SUBROUTINE SMALLTICK(XX1,XX2,YY1,YY2,XLOW,XHIGH,SPAN,
     C                     XTLEN,TYPE)
C
C PRODUCE SMALL TICK MARKS ALONG INDEPENDENT AXIS (TIME)
C AT USER DEFINED INTERVALS IN DAYS.
C
      INTEGER*4 TYPE
      REAL*8 XLOW,XHIGH,XVAL,VAL,YR1,YR2,FRAC1,FRAC2
C
      Y1 = YY1 + XTLEN/2.
      Y2 = YY2 - XTLEN/2.
      IF (TYPE.GT.-5) THEN
       XDIFF = REAL((XHIGH-XLOW)*365.25D0)
      ELSE
       XDIFF = REAL(XHIGH-XLOW)
      END IF
      YDIFF = YY2 - YY1
      ISTOP = NINT(XDIFF/SPAN) + 1
      XVAL = XLOW
      DO 50 I=1,ISTOP
       XLIN = XX1 + (XX2-XX1) * REAL((XVAL-XLOW)/(XHIGH-XLOW))
       IF (XVAL.GT.XHIGH) GOTO 50
       CALL MGOLINE(XLIN,YY1,XLIN,Y1)
       CALL MGOLINE(XLIN,YY2,XLIN,Y2)
       IF (TYPE.GT.-5) THEN
        YR1 = DBLE(INT(XVAL))
        VAL = XVAL + DBLE(SPAN) / DBLE(DAYS(XVAL))
        IF (INT(VAL).GT.INT(XVAL)) THEN
         YR2 = YR1 + 1.D0
         FRAC1 = (YR2 - XVAL) * DBLE(DAYS(XVAL))
         FRAC2 = DBLE(SPAN) - FRAC1 
         VAL = YR2 + FRAC2 / DBLE(DAYS(YR2))
        ENDIF
        XVAL = VAL
       ELSE
        XVAL = XVAL + SPAN
       END IF
   50 CONTINUE
C
      RETURN
      END
      SUBROUTINE LRGTICK(VV1,VV2,WW1,WW2,VLOW,VHIGH,
     C                   SPAN,FORM,TLEN,LAB,ISCALE)
      CHARACTER*1 TCHAR
      CHARACTER*5 EXPONENT
      CHARACTER*10 FORM,NCHAR
      REAL*4 NUMH
C
C DRAW LARGE TICK MARKS AND LABEL
C
      IF (VV1.LT.0) THEN
       VV1 = ABS(VV1)
       IAX = 2
      ELSE 
       IAX = 1
      END IF
      DIFF = VHIGH - VLOW
      ISTOP = NINT(DIFF/SPAN) + 1
      V0 = VV1 - TLEN/2.
      V1 = VV1 + TLEN
      V2 = VV2 - TLEN
      VAL = VLOW
      DO I=1,ISTOP 
       VLIN = WW1 + (WW2-WW1)*(VAL-VLOW)/DIFF
       ERR = (VLIN-WW2)/(WW2-WW1)
       IF (ERR.GT.0.0001) GOTO 100
       IF (IAX.EQ.1) THEN
C
C DRAW INDEPENDENT AXIS LARGE TICKS
C
        CALL MGOLINE(VLIN,VV1,VLIN,V1)
        CALL MGOLINE(VLIN,VV2,VLIN,V2)
       ELSE IF (IAX.EQ.2) THEN
C
C DRAW DEPENDENT AXIS LARGE TICKS
C
        CALL MGOLINE(VV1,VLIN,V1,VLIN)
        CALL MGOLINE(VV2,VLIN,V2,VLIN)
       END IF
C
       IF ( LAB.EQ.4 ) GOTO 60
C
C LINEAR SCALING
C
       IF (ISCALE.EQ.0) THEN
        TCHAR = FORM(1:1)
        IF (TCHAR.EQ.'I') THEN
         WRITE(NCHAR,'('//FORM//')') INT(VAL)
        ELSE
         WRITE(NCHAR,'('//FORM//')') VAL
        END IF
       END IF
C
C LOGARITHMIC SCALING
C
       IF (ISCALE.EQ.1) THEN
        WRITE(EXPONENT,'(I5)') INT(VAL)
        WRITE(NCHAR,'(A10)') '10\\U'//EXPONENT
       END IF
C
C SHIFT LEFT CHARACTER STRING
C
       NUM = 0
       DO 55 ILET=1,10
        TCHAR = NCHAR(ILET:ILET)
        LVAL = ICHAR(TCHAR)
        IF ( LVAL.LT.32 .OR. LVAL.EQ.127 ) GOTO 55
        IF ( LVAL.NE.32 ) THEN
         NUM = NUM + 1
         NCHAR(NUM:NUM) = TCHAR
         IF (NUM.NE.ILET) NCHAR(ILET:ILET) = CHAR(32)
        END IF
   55  CONTINUE
C
       IF ( LAB.EQ.1.AND.I.EQ.1 ) GOTO 60
       IF ( LAB.EQ.2.AND.I.EQ.ISTOP ) GOTO 60
       IF ( LAB.EQ.3.AND.(I.EQ.1.OR.I.EQ.ISTOP) ) GOTO 60
       IF (IAX.EQ.1) THEN
C        TYPE*,'X LABEL *'//NCHAR//'*'
        CALL MGOGRELOCATE(VLIN,V0)
        CALL MGOPUTLABEL(NUM,NCHAR,2)
       END IF
       IF (IAX.EQ.2) THEN
C        TYPE*,'Y LABEL *'//NCHAR//'*'
        CALL MGOGRELOCATE(V0,VLIN)
        CALL MGOPUTLABEL(NUM,NCHAR,4)
       END IF
   60  VAL = VAL + SPAN
      END DO
  100 CONTINUE
      RETURN 
      END
      SUBROUTINE SMLTICK(VV1,VV2,WW1,WW2,VLOW,VHIGH,SPAN,TLEN)
C
C DRAW SMALL TICK MARKS
C
      IF (VV1.LT.0) THEN
       VV1 = ABS(VV1)
       IAX = 2
      ELSE 
       IAX = 1
      END IF
      DIFF = VHIGH - VLOW
      ISTOP = NINT(DIFF/SPAN) + 1
      V1 = VV1 + TLEN/2.
      V2 = VV2 - TLEN/2.
      VAL = VLOW
      DO I=1,ISTOP
       VAL = VAL + SPAN
       VLIN = WW1 + (WW2-WW1)*(VAL-VLOW)/DIFF
       IF (VLIN.GT.WW2) GOTO 45
       IF (IAX.EQ.1) THEN
C
C DRAW INDEPENDENT AXIS SMALL TICKS
C
        CALL MGOLINE(VLIN,VV1,VLIN,V1)
        CALL MGOLINE(VLIN,VV2,VLIN,V2)
       ELSE IF (IAX.EQ.2) THEN
C
C DRAW DEPENDENT AXIS SMALL TICKS
C
        CALL MGOLINE(VV1,VLIN,V1,VLIN)
        CALL MGOLINE(VV2,VLIN,V2,VLIN)
       END IF
C
      END DO
   45 CONTINUE
      RETURN
      END
      SUBROUTINE LOGTICK(VV1,VV2,WW1,WW2,VLOW,VHIGH,SPAN,TLEN)
C
C DRAW SMALL TICK MARKS
C
      IF (VV1.LT.0) THEN
       VV1 = ABS(VV1)
       IAX = 2
      ELSE 
       IAX = 1
      END IF
      DIFF = VHIGH - VLOW
      ISTOP = NINT(DIFF/SPAN) + 1
      V1 = VV1 + TLEN/2.
      V2 = VV2 - TLEN/2.
      DO DECADE = VLOW,VHIGH,1.0
       CURDEC = DECADE
       DO I = 1,10
        VAL = ALOG10(REAL(I)*10.0**CURDEC)
        VLIN = WW1 + (WW2-WW1)*(VAL-VLOW)/DIFF
        IF (VLIN.GT.WW2) GOTO 45
        IF (IAX.EQ.1) THEN
C
C DRAW INDEPENDENT AXIS SMALL TICKS
C
         CALL MGOLINE(VLIN,VV1,VLIN,V1)
         CALL MGOLINE(VLIN,VV2,VLIN,V2)
        ELSE IF (IAX.EQ.2) THEN
C
C DRAW DEPENDENT AXIS SMALL TICKS
C
         CALL MGOLINE(VV1,VLIN,V1,VLIN)
         CALL MGOLINE(VV2,VLIN,V2,VLIN)
        END IF
C
       END DO
      END DO
   45 CONTINUE
      RETURN
      END
      SUBROUTINE CALENDAR(TIME,IYEAR,IDAY,IHOUR)
C 
C CONVERT DOUBLE PRECISION DECIMAL YEAR INTO INTEGER YEAR, DAY AND HOUR 
C

      REAL*8 TIME,DYEAR,LEAP
C
      IYEAR = INT(TIME)
      DYEAR = DBLE(IYEAR)
      LEAP = DBLE(DAYS(TIME))
      DDAY = REAL((TIME-DYEAR)*LEAP) + 1.
      IDAY = INT(DDAY)
      DHOUR = (DDAY-REAL(IDAY))*24.
      IHOUR = INT(DHOUR)
      DMIN = DHOUR - REAL(IHOUR)
      IF (DMIN.GT.0.8) IHOUR = IHOUR + 1
C
      IF (IHOUR.EQ.24) THEN
       IDAY = IDAY + 1
       IHOUR = 0
      END IF
C
      ILEAP = INT(LEAP)
      IF (IDAY.GT.ILEAP) THEN
       IYEAR = IYEAR + 1
       IDAY = IDAY - ILEAP
      END IF
C
      RETURN
      END
      REAL*8 FUNCTION CONVERT(XVAL)
      CHARACTER*10 VALUE
      REAL*8 DEC
C
C GET START AND STOP TIMES
C CONVERT CONCATENATED DATE INTO DECIMAL YEAR
C
      IVAL = INT(XVAL)
      WRITE(VALUE,'(I7)') IVAL
      READ(VALUE,'(I2,I3.3,I2.2)') IY,ID,IH
      CONVERT = DEC(IY,ID,IH)
      RETURN
      END
      REAL*8 FUNCTION DEC(NYR,NDY,NHR)
C
C CONVERT CALENDAR TIME INTO DECIMAL YEAR
C
      REAL*8 YEAR,DAY,HOUR
      YEAR = DFLOAT(NYR)
      DAY = DFLOAT(NDY)
      HOUR= DFLOAT(NHR)
      DEC = YEAR + (DAY-1.D0 + HOUR/24.D0) / DBLE(DAYS(YEAR))
C
      RETURN
      END
      REAL*4 FUNCTION DAYS(TIME)
C
C COMPUTE LENGTH OF YEAR IN DAYS
C
      REAL*8 TIME
      NYR = INT(TIME)
      DAYS = 365.
      IF (MOD(NYR,4).EQ.0) DAYS = 366.
C
      RETURN
      END
      SUBROUTINE MCONNECT(FILL,XVAL,YVAL,PTS,DELTA)
C
C ROUTINE TO CONNECT ADJACENT ARRAY POINTS WITH INDEPENDENT DELTAS
C LESS THAN <DELTA> WHILE REMOVING FILL VALUES.
C
      INTEGER*4 PTS
      REAL*4 X(2),Y(2),XVAL(1),YVAL(1)
C
C FILL  FILL DATA VALUE
C XVAL  INDEPENDENT AXIS ARRAY
C YVAL  DEPENDENT AXIS ARRAY
C PTS   NUMBER OF ELEMENTS IN PARALLEL ARRAYS XVAL & VAL
C DELTA MAXIMUM GAP BETWEEN TWO ADJACENT POINTS TO BE CONNECTED 
C
C      WRITE(6,*)
C      WRITE(6,*) '***ENTER PLOTTING SUBROUTINE***'
      IPT = 0
      N = PTS - 1
      DO 150 I=1,N
       IF ( YVAL(I).NE.FILL .AND. YVAL(I+1).NE.FILL .AND.
     .      (XVAL(I+1) - XVAL(I)).LT.DELTA ) THEN
        X(1) = XVAL(I)
        X(2) = XVAL(I+1)
        Y(1) = YVAL(I)
        Y(2) = YVAL(I+1)
        CALL MGOCONNECT(X,Y,2)
        IPT = IPT + 1
       ELSE IF ( YVAL(I).NE.FILL ) THEN
        CALL MGORELOCATE(XVAL(I),YVAL(I))
        CALL MGOPOINT(10,3)
        IPT = IPT + 1
        IF ( I.EQ.N.AND.YVAL(I+1).NE.FILL ) THEN
         CALL MGORELOCATE(XVAL(I+1),YVAL(I+1))
         CALL MGOPOINT(10,3)
         IPT = IPT + 1
        END IF
       END IF
  150 CONTINUE
C      WRITE(6,*) 'NUMBER OF INTERVALS/POINTS : ',IPT
C      WRITE(6,*) '***EXIT PLOTTING SUBROUTINE***'
C
      RETURN
      END
      SUBROUTINE GETDATE(SYSDATE)
C
C CALL DEC FORTRAN SUBROUTINE TO GET SYSTEM DATE 
C
      CHARACTER SYSDATE*9
C
      CALL DATE(SYSDATE)
C
      RETURN 
      END
      SUBROUTINE GETTIME(SYSTIME)
C
C CALL DEC FORTRAN SUBROUTINE TO GET SYSTEM TIME
C
      CHARACTER SYSTIME*8
C
      CALL TIME(SYSTIME)
C
      RETURN
      END