      PROGRAM TEST
C
C     PURPOSE: ILLUSTRATE USAGE OF SUBROUTINES HEIGHTDRIVER, HEIGHT, AND
C              ASSOCIATED SUBROUTINES
C              THIS PROGRAM DUPLICATES THE FIRST 10 ENTRIES
C              (TIME,LAT,LON,HEIGHT)
C              IN THE GTRACK FILE (GRNDTK). THIS COULD BE USED
C              AS A FIRST TEST TO MAKE SURE THE CODE WORKS
C
C     GERARD L.H. KRUIZINGA                   10/20/93
C
C     MODIFIED BY: GERARD L.H. KRUIZINGA      05/08/94
C                  MODIFICATION TO TAKE INTO ACCOUNT PRECISION
C                  OTHER THAN CRAY 64-BIT MACHINES
C     MODIFIED BY: MICHAEL J. GABOR  3/20/95
C                  MODIFICATION TO TAKE BURN TIMES INTO ACCOUNT
C

      IMPLICIT NONE

      DOUBLE PRECISION EPOCHORB,EPOCHREC
      DOUBLE PRECISION DEPOCH,DEPOCH1,EPOCH,TORB(4)
      DOUBLE PRECISION TIME,DELTT,TBOFF(1000),TBON(1000)
      DOUBLE PRECISION Y(4,3),HT,QLAT,QLON

      DOUBLE PRECISION SECONDS,AE,RFLAT,DELT

      INTEGER          ISATID,IDATE(5)
      INTEGER          INITFLG,IUNITG,IUNITO,I
      INTEGER          IFLAG,IUNITM,NUMBRN

      CHARACTER*10     JOBTIM,JOBDAT,NAMESAT,VERSION

      LOGICAL MFLAG

      COMMON /HTSTUFF/ EPOCHORB,DEPOCH,DEPOCH1,TORB,Y,INITFLG
      COMMON /GTRINFO/ DELT,ISATID,IDATE,AE,RFLAT,EPOCH,SECONDS,
     .                 VERSION,JOBDAT,JOBTIM,NAMESAT
      DATA MFLAG/.TRUE./

      DELTT=60.D0
      IUNITG  = 11
      IUNITO  = 12
      IUNITM  = 13

C  INITIALIZE NUMBRN BEFORE THE FIRST CALL OF HEIGHTDRIVER SUBROUTINE.
C	SUBROUTINE UPDATES IT TO ACTUAL VALUE.  MAKE SURE THAT INITIALLY
C	IT IS LARGE

      NUMBRN=1000

	OPEN(IUNITO,FILE='OUTPUT')
	OPEN(IUNITG,FILE='GRNDTK',STATUS='OLD')
	OPEN(IUNITM,FILE='MANEUVERS',STATUS='OLD')
	REWIND(IUNITM)

      CALL CHECK_HEADER(IUNITG)

      REWIND (IUNITG)
C
C     SOME COMMENTS ON CALLING SUBROUTINE HEIGHTDRIVER
C
C     VARIABLE "TIME" IS IN SECONDS PAST EPOCHREC (IN MJD)
C     VARIABLE "HT"   IS IN METERS
C     VARIABLES "QLAT" AND "QLON" ARE IN DEGREES (QLAT IS GEODETIC LATITUDE)
C     VARIABLE "DELTT" SPECIFIES REGION OF PERMITTED TIME AROUND BURN IN
C		SECONDS.  TOTAL REGION IS DELTT ON EITHER SIDE OF BURN.
C     ARRAY  "TBON" AND "TBOFF" SPECIFY THE START AND STOP TIMES OF THE
C		BURNS AND ARE AVAILABLE FROM THE MANEUVERS TABLE
C

      EPOCHREC= EPOCH

      DO 10 I = 110530,110570  
C     DO 10 I = 1,20  
          TIME= DFLOAT(I-1)*10.0D0
          CALL HEIGHTDRIVER (IUNITM,DELTT,MFLAG,NUMBRN,TBON,TBOFF,
     &               TIME,EPOCHREC,IUNITG,HT,QLAT,QLON,IFLAG,IUNITO)
          WRITE(IUNITO,1000) TIME,QLAT,QLON,HT
10    CONTINUE

       CLOSE(IUNITG)
       CLOSE(IUNITO)
       CLOSE(IUNITM)
1000	FORMAT (F15.6,2F15.8,F15.4)
      STOP
      END

        SUBROUTINE HEIGHTDRIVER (IUNITM,DELTT,MFLAG,NUMBRN,TBON,TBOFF,
     &          TIME,EPOCHREC,IUNITG,HT,QLAT,QLON,IFLAG,IUNITO)
************************************************************************
*  MODIFIED BY: MICHAEL J. GABOR  3/20/95
*	        MODIFICATION TO TAKE BURN TIMES INTO ACCOUNT
*  THIS IS A MODIFICATION TO THE HEIGHT SUBROUTINE (DETAILS UNDER
*     HEIGHT) THAT ALLOWS FOR DATA TO BE IGNORED THAT IS WITHIN
*     A SPECIFIED NUMBER OF SECONDS OF A BURN.  DATA IS ARBITRARILY
*     REMOVED DURING THE ACTUAL BURN.  TO INCLUDE DATA THAT IS DURING
*     A BURN, SPECIFY 0.D0 SECONDS FOR DELTT.
*  VARIABLES SPECIFICALLY USED IN HEIGHTDRIVER:
*     ARRAY  "TBON" AND "TBOFF" SPECIFY THE START AND STOP TIMES OF THE
*               BURNS AND ARE AVAILABLE FROM THE MANEUVERS TABLE
*     VARIABLE "MFLAG" DETERMINES IF THIS IS THE FIRS PASS THROUGH THE
*		LIST OF BURNS
*     VARIABLE "DELTT" SPECIFIES REGION OF PERMITTED TIME AROUND BURN IN
*               SECONDS.  TOTAL REGION IS DELTT ON EITHER SIDE OF BURN.
************************************************************************

	IMPLICIT NONE
	DOUBLE PRECISION DELTT,TBON(1000),TBOFF(1000),TIME,EPOCHREC,HT
	DOUBLE PRECISION QLAT,QLON
	INTEGER IUNITM,IUNITG,NUMBRN,I,IFLAG,IUNITO
	LOGICAL MFLAG,STATUS
C       DATA STATUS/.TRUE./
	SAVE

        STATUS = .TRUE.
        QLAT   = -999.0D0
        QLON   = -999.0D0
        HT     = -999.0D0

	DO 100 I=1,NUMBRN
		IF (MFLAG.EQ..TRUE.) THEN
		    READ(IUNITM,5,END=200) TBON(I),TBOFF(I)
		    TBON(I)=(TBON(I)-2400000.5D0-EPOCHREC)*86400.D0
	   	    TBOFF(I)=(TBOFF(I)-2400000.5D0-EPOCHREC)*86400.D0
		    NUMBRN=I
		END IF
	
		IF (DABS(TIME-TBOFF(I)).LT.DELTT.OR.DABS(TIME-TBON(I))
     &			.LT.DELTT) THEN
			STATUS=.FALSE.
		END IF
	
		IF (TIME.GT.TBON(I).AND.TIME.LT.TBOFF(I).AND.DELTT.GT.0.D0)
     &			THEN
			STATUS=.FALSE.
		END IF

100	CONTINUE
200	MFLAG=.FALSE.

	IF (STATUS.EQ..TRUE.) THEN
		CALL HEIGHT (TIME,EPOCHREC,IUNITG,HT,QLAT,QLON,IFLAG)
	ELSE
	END IF

5	FORMAT(2X,F18.10,4X,F18.10)
1000	FORMAT (F15.6,2F15.8,F15.4)
	END


      SUBROUTINE HEIGHT (TREC,EPOCHREC,UNIT,HT,LAT,LONG,IFLAG)
************************************************************************
*  CODED BY: DON CHAMBERS -- 2/4/92 -- VERSION 1.3
* 
*  MODIFIED BY: GERARD L.H. KRUIZINGA  4/1/92
*               INCLUSION OF CORRECT SPLINE INTERPOLATION ACROSS
*               THE PRIME (GREENWHICH MERIDIAN)
*  MODIFIED BY: GERARD L.H. KRUIZINGA  5/8/94
*               MODIFICATIONS TO TAKE INTO ACCOUNT PRECISION
*               DIFFERENCES BETWEEN CRAY AND SUN WORKSTATIONS 
* 
*  THIS SUBROUTINE READS THROUGH A UTOPIA GTRACK FILE CONTAINING HEIGHTS
*  ABOVE THE REFERENCE ELLIPSOID SPACED AT 10 SECOND INTERVALS FROM THE
*  EPOCH OF THE CYCLE.  
*
*  THIS SUBROUTINE INTERPOLATES (USING A CUBIC SPLINE) TO FIND THE ORBIT
*  HEIGHT AT TREC SO THAT THE ORBIT HEIGHT CAN BE COMPARED TO THE GDR
*  HEIGHT.
*  
*
*     VARIABLES:
*  INPUT--
*     TREC - TIME OF RECORD (SECONDS PAST EPOCHREC)
*     EPOCHREC - EPOCH OF RECORD (IN MJD)
*     UNIT - UNIT NUMBER FOR ORBIT HEIGHTS FILE
*  OUTPUT--
*     HT - ORBITAL HEIGHT ABOVE REFERENCE ELLIPSOID (M)
*     LAT - GEODETIC LATITUDE OF SUB-SATELLITE POINT (IN DEGREES)
*     LONG - GEOCENTRIC LONGITUDE OF SUB-SATELLITE POINT (IN DEGREES)
*     IFLAG - ERROR FLAG 
*             IFLAG = 0 : SUCCESFUL INTERPOLATION
*             IFLAG =-1 : REQUESTED TIME BEFORE EPOCH
*             IFLAG = 1 : REQUESTED TIME AFTER LAST RECORD IN FILE 
*  OTHERS--
*     EPOCHORB - EPOCH OF ORBIT CYCLE (IN JD)
*     TORB - TIME OF ORBIT PAST EPOCHORB (IN SECONDS)
*     HTERR - ERROR MESSAGES
*     
*  MODIFICATION--
*     HEIGHT SUBROUTINE CAN READ GTRACK FILES WITH AND WITHOUT   
*     GTRACK HEADERS. SUBROUTINE WILL CHECK IF CORRECT TIME SPACING
*     ON APPEARS ON GTRACK HEADER
*
*  MODIFICATION
*     INTERPOLATION WILL TAKE INTO ACCOUNT LEAP SECONDS 
* 
****************************************************************************

      IMPLICIT NONE

      DOUBLE PRECISION TREC,EPOCHORB,EPOCHREC
      DOUBLE PRECISION DEPOCH,DEPOCH1,TORB(4),TOLD,EPOCH,DTNEW
      DOUBLE PRECISION LAT,LONG
      DOUBLE PRECISION SECONDS,AE,RFLAT,DELT
      DOUBLE PRECISION Y(4,3),HT
      DOUBLE PRECISION DT,DY2,DY3,DY4,A,B,C,VAL

      INTEGER          UNIT
      INTEGER          ISATID,IDATE(5)
      INTEGER          INITFLG,I,J,IFLAG

      CHARACTER*80     RECORD,HTERR
      CHARACTER*10     JOBTIM,JOBDAT,NAMESAT,VERSION
      CHARACTER*5      TEXT

      LOGICAL          FOUND,FLAG

      COMMON /HTSTUFF/ EPOCHORB,DEPOCH,DEPOCH1,TORB,Y,INITFLG
      COMMON /GTRINFO/ DELT,ISATID,IDATE,AE,RFLAT,EPOCH,SECONDS,
     .                 VERSION,JOBDAT,JOBTIM,NAMESAT


      DATA  NAMESAT  /'NOTDET'/
      DATA  TOLD     /999999999999999.0/

      FOUND = .FALSE.
      HTERR = 'ERROR IN READING ORBIT HEIGHT FILE'
      IFLAG = 0

      DO WHILE (.NOT. FOUND)
*     
*     READ IN FIRST 4 POINTS THE FIRST TIME THE SUBROUTINE IS CALLED
*     
         IF (INITFLG .EQ. 0) THEN
            TOLD = 99999999999999.0
            CALL CHECK_HEADER(UNIT)
            IF (NAMESAT .EQ. 'NOTDET') THEN
               READ (UNIT,1000,END=99,ERR=199) RECORD
               READ(RECORD,1200,ERR=299) TEXT,EPOCHORB
               DEPOCH  = (EPOCHORB - EPOCHREC)*86400.0D0
               DEPOCH1 = DEPOCH
            ELSE
               IF (DELT .NE. 10.0) THEN 
                  WRITE(*,*) ' INCORRECT TIME INTERVAL FOR HEIGHT'
                  WRITE(*,*) ' INTERVAL : ',DELT,' SHOULD BE 10.0 SEC'
                  STOP
               ENDIF
               DEPOCH  = (EPOCH - EPOCHREC)*86400.0D0
               DEPOCH1 = DEPOCH
            ENDIF           
            DO 10 J =1,4
 9             READ (UNIT,1000,END=99,ERR=199) RECORD
               READ (RECORD,1100)TORB(J),Y(J,1),Y(J,2),Y(J,3)
               TORB(J) = TORB(J) + DEPOCH
               IF (TOLD .EQ. TORB(J)) GOTO 9
               TOLD = TORB(J)
 10         CONTINUE
            INITFLG = 1
            WRITE(*,*) ' '         
            WRITE(*,*) 'INITIALIZATION INFORMATION FROM HEIGHT'
            WRITE(*,*) ' ' 
            IF (NAMESAT .EQ. 'NOTDET') THEN        
              WRITE(*,*) 'EPOCH READ FROM GTRACK FILE :',EPOCHORB
            ELSE
              WRITE(*,*) 'EPOCH READ FROM GTRACK FILE :',EPOCH
            ENDIF             
            WRITE(*,*) 'EPOCH FOR TIME TAG          :',EPOCHREC
            WRITE(*,*) 'RELATIVE EPOCH              :',DEPOCH
            WRITE(*,*) ' '             
         END IF


      IF (TREC.LT.DEPOCH1) THEN
         HTERR = 'RECORD TIME IS BEFORE ORBIT TIME'
         IFLAG=-1
         GO TO 199
      END IF
*
* SEE IF T LIES IN THE FOUR POINTS; IF SO, INTERPOLATE TO FIND
* HT, LAT, AND LONG; IF NOT, SEARCH THROUGH NEXT FOUR POINTS
*
      IF (TREC.GE.TORB(1) .AND. TREC.LT.TORB(4)) THEN
         DO 30 I =1,3
            CALL CHECK_LONGITUDE
            DY2 = Y(2,I) - Y(1,I)
            DY3 = Y(3,I) - Y(1,I)
            DY4 = Y(4,I) - Y(1,I)
            DT = TREC - TORB(1)
            CALL LEAPSEC(TREC,TORB,EPOCHREC,DY2,DY3,DY4,
     &                   A,B,C,FLAG,DTNEW)
            IF (FLAG) THEN
              DT=DTNEW
              GOTO 29
            ENDIF
            A = (DY4 + 3.0*DY2 - 3.0*DY3)/6000.0
            B = (4.0*DY3 - 5.0*DY2 - DY4)/200.0
            C = (DY4 - 4.5*DY3 + 9.0*DY2)/30.0
 29         VAL = A*DT*DT*DT + B*DT*DT + C*DT + Y(1,I)
            IF (I.EQ.1) LAT = VAL
            IF (I.EQ.2) LONG = VAL
            IF (I.EQ.2 .AND. VAL .GT. 360.0) LONG = VAL - 360.0
            IF (I.EQ.3) HT = VAL
 30      CONTINUE
         IF (Y(4,2) .GT. 360.0) Y(4,2)=Y(4,2)-360.0
         FOUND = .TRUE.
      ELSE
         TORB(1) = TORB(4)
         Y(1,1) = Y(4,1)
         Y(1,2) = Y(4,2)
         Y(1,3) = Y(4,3)
         DO 35 J = 2,4
 721     READ (UNIT,1000,END=99,ERR=199) RECORD
         IF (RECORD(1:5) .EQ. 'EPOCH') THEN
            READ(RECORD,1200) TEXT,EPOCHORB
            DEPOCH = (EPOCHORB - EPOCHREC)*86400.0
*
* SKIP FIRST RECORD AFTER EACH EPOCH CARDS PAST THE FIRST ONE
*
            READ (UNIT,1000,END=99,ERR=199) RECORD   
            READ (UNIT,1000,END=99,ERR=199) RECORD   
         END IF
         READ (RECORD,1100)TORB(J),Y(J,1),Y(J,2),Y(J,3)
         TORB(J) = TORB(J) + DEPOCH
         IF (TOLD .EQ. TORB(J)) GOTO 721
         TOLD = TORB(J)
 35      CONTINUE
      END IF
      END DO
      IF (.NOT. FOUND) THEN
         HTERR = 'RECORD TIME IS BIGGER THAN LARGEST ORBIT TIME'
         IFLAG=1
         GO TO 199
      END IF
      RETURN
C 99   PRINT *,'END OF ORBIT HEIGHT FILE REACHED'
C      PRINT *,'LAST REQUESTED TIME :',TREC
 99   CONTINUE
      IFLAG=1
      RETURN
C 199  PRINT *,'ERROR OCCURRED IN SUBROUTINE HEIGHT'
C      PRINT *,HTERR
 199  CONTINUE
      RETURN
 299  PRINT *,'ERROR READING FIRST LINE OF GTRACK FILE; NO EPOCH HEADER'
      STOP
 1000 FORMAT (A80)
 1100 FORMAT (F15.6,2F15.8,F15.4)
 1150 FORMAT (F15.7,2F15.8,F15.4)
 1200 FORMAT (A5,1X,F25.15)
      END
      SUBROUTINE CHECK_LONGITUDE
*     
*     SUBROUTINE TO CHECK IF ALL 4 LONGITUDES DO NOT CONTAIN A 
*     DISCONTINUITY CAUSED BY CROSSING THE PRIME MERIDIAN
*     CODED BY : GERARD L.H. KRUIZINGA     4/2/92

      IMPLICIT NONE

      DOUBLE PRECISION EPOCHORB,DEPOCH,DEPOCH1,TORB(4)
      DOUBLE PRECISION Y(4,3)

      INTEGER          INITFLG,I

      COMMON /HTSTUFF/ EPOCHORB,DEPOCH,DEPOCH1,TORB(4),Y(4,3),INITFLG


      IF (ABS(Y(1,2)-Y(4,2)).GT. 100.0) THEN
        DO 10 I=1,4
          IF ( Y(I,2) .LT. 100) Y(I,2)=Y(I,2)+360.0
 10     CONTINUE
      ENDIF
      RETURN
      END
C**********************************************************************
      SUBROUTINE KALDAY (ITIME, STIME, UTC)
C
C...PURPOSE:  TO CONVERT BETWEEN CALENDAR DATE AND JULIAN DATE
C
C...PROGRAMMED BY JAMES D MCMILLAN  -  UNIV. OF TEXAS  -  7/23/1973
C   JULDAY ADDED BY BRIAN CUTHBERTSON  -  UNIV. OF TEXAS  -  9/29/1979
C
C...KALDAY INPUT  (JULDAY OUTPUT):
C      UTC       THE DOUBLE PRECISION JULIAN DATE 
C
C...KALDAY OUTPUT  (JULDAY INPUT):
C      ITIME   INTEGER ARRAY WITH MONTH, DAY, YEAR, HOUR AND MINUTE
C      STIME   SINGLE PRECISION SECONDS 
C
      DIMENSION ITIME(5)
      DOUBLE PRECISION UTC
C
C...COMPUTE THE HOUR, MINUTE, AND SECOND.
      JD=UTC
      XTIME=AMOD(SNGL(UTC-JD)*24.+12.,24.)
      ITIME(4)=XTIME
      XTIME=(XTIME-ITIME(4))*60.0
      ITIME(5)=XTIME
      STIME=(XTIME-ITIME(5))*60.0
C
C...ROUND 60 SECONDS TO NEXT MINUTE
      IF (STIME.LT.59.99999999999) GO TO 10
      STIME=0.0
      ITIME(5)=ITIME(5)+1
      IF (ITIME(5).NE.60) GO TO 10
      ITIME(5)=0
      ITIME(4)=ITIME(4)+1
      IF (ITIME(4).NE.24) GO TO 10
      ITIME(4)=0
C
C...COMPUTE THE YEAR, MONTH, AND DAY.
   10 IF (ITIME(4).LT.12) JD=JD+1
      LX=JD+68569
      NX=4*LX/146097
      LX=LX-(146097*NX+3)/4
      ITIME(3)=4000*(LX+1)/1461001
      LX=LX-1461*ITIME(3)/4+31
      ITIME(1)=80*LX/2447
      ITIME(2)=LX-2447*ITIME(1)/80
      LX=ITIME(1)/11
      ITIME(1)=ITIME(1)+2-12*LX
      ITIME(3)=100*(NX-49)+ITIME(3)+LX
      RETURN
C
      ENTRY JULDAY (ITIME, STIME, UTC)
C
      JD=ITIME(2)-32075+1461*(ITIME(3)+4800+(ITIME(1)-14)/12)/4+367*(ITI
     1ME(1)-2-(ITIME(1)-14)/12*12)/12-3*((ITIME(3)+4900+(ITIME(1)-14)/12
     2)/100)/4
      UTC=DBLE(FLOAT(JD))-0.5D+00+DBLE(FLOAT(ITIME(4)))/24.0D+00+DBLE(FL
     1OAT(ITIME(5)))/1440.0D+00+DBLE(STIME)/86400.0D+00
      RETURN
C
C.....END KALDAY/JULDAY
      END 
C**********************************************************************
      SUBROUTINE GET_UTC(KYR,KDAY,KSEC,KMSEC,TGDR)
C
C     THIS SUBROUTINE COMPUTES GEOSAT GDR TIME (SEC PAST 01/01/1985)
C
C     CODED BY GERARD L.H. KRUIZINGA             JULY 199OA2
C     CENTER FOR SPACE RESEARCH           UNIVERSITY OF TEXAS AT AUSTIN
C
C     INPUT    : KYR    YEAR
C                KDAY   DAY IN THE YEAR
C                KSEC   SECONDS PAST MIDNIGHT OF KDAY
C                KMSEC  MICRO SECONDS
C     OUTPUT   : TGDR   SECONDS PAST 01/01/1985
C
      DOUBLE PRECISION TGDR,UTC,JAN1,GEOSAT_EPOCH
      DIMENSION ITIME(5)
      DATA ITEST /0/
      DATA GEOSAT_EPOCH /0.0D0/

      IF (ITEST .EQ. 0) THEN
        ITIME(1) =  1
        ITIME(2) =  1
        ITIME(3) =  1985
        ITIME(4) =  0
        ITIME(5) =  0
        STIME    =  0.0
        CALL JULDAY (ITIME,STIME,UTC)
        GEOSAT_EPOCH = UTC
        ITEST = 1
      ENDIF

      ITIME(1) =  1
      ITIME(2) =  1
      ITIME(3) =  1900 + KYR
      ITIME(4) =  0
      ITIME(5) =  0
      STIME    =  0.0

      CALL JULDAY (ITIME, STIME, JAN1)

      TGDR = (JAN1 - GEOSAT_EPOCH)*86400.0 + DBLE(KDAY-1)*86400.0
     &       + DBLE(KSEC) + FLOAT(KMSEC)/1E6

      RETURN
      END
      SUBROUTINE CHECK_HEADER(IUNIT)
C     
C     PURPOSE: CHECK FOR THE EXISTANCE OF AN GTRACK HEADER 
C              AND IF SO READ THE INFO
C     
C     CODED BY: GERARD L.H. KRUIZINGA       05/07/93
C     CENTER FOR SPACE RESEARCH       UNIVERSITY OF TEXAS AT AUSTIN
C     
      IMPLICIT NONE

      DOUBLE PRECISION  EPOCH

      DOUBLE PRECISION  SECONDS,AE,RFLAT,DELT

      INTEGER           ISATID,IDATE(5),IUNIT

      CHARACTER*10      JOBTIM,JOBDAT,NAMESAT,VERSION
      CHARACTER*60      CARD

      COMMON /GTRINFO/ DELT,ISATID,IDATE,AE,RFLAT,EPOCH,SECONDS,
     .                 VERSION,JOBDAT,JOBTIM,NAMESAT

      READ(IUNIT,390) CARD
      READ(IUNIT,390) CARD
      READ(IUNIT,390) CARD
      READ(IUNIT,390) CARD
      PRINT*,'CARD = ',CARD

      IF (CARD(1:7) .EQ. 'VERSION') THEN
        IF (CARD(1:8) .EQ. 'VERSION:') THEN
          READ(CARD,400) VERSION,DELT
          READ(IUNIT,410) JOBDAT,JOBTIM,ISATID,NAMESAT
          READ(IUNIT,390) CARD
          READ(IUNIT,420) AE
          READ(IUNIT,430) RFLAT
          READ(IUNIT,390) CARD
          READ(IUNIT,440) EPOCH,IDATE,SECONDS
          READ(IUNIT,390) CARD
          READ(IUNIT,190) CARD
          READ(IUNIT,390) CARD
        ELSE
          READ(CARD,200) VERSION,DELT
          READ(IUNIT,210) JOBDAT,JOBTIM,ISATID,NAMESAT
          READ(IUNIT,190) CARD
          READ(IUNIT,220) AE
          READ(IUNIT,230) RFLAT
          READ(IUNIT,190) CARD
          READ(IUNIT,240) EPOCH,IDATE,SECONDS
          READ(IUNIT,190) CARD
          READ(IUNIT,190) CARD
          READ(IUNIT,190) CARD
        ENDIF
      ELSE
        REWIND(IUNIT)
      ENDIF
      WRITE(*,*) DELT,ISATID,AE,RFLAT,EPOCH,IDATE,SECONDS,
     .                 VERSION,JOBDAT,JOBTIM,NAMESAT

 190  FORMAT (A79)
 200  FORMAT (10X,A10,30X,F9.1)
 210  FORMAT (9X,2A10,23X,I7,1X,A10)
 220  FORMAT (49X,F10.1)
 230  FORMAT (49X,F10.5)
 240  FORMAT (8X,F20.12,12X,I3,1X,I2,1X,I2,I4,1X,I2,1X,F10.7)
 390  FORMAT (A60)
 400  FORMAT (9X,A10,22X,F9.1)
 410  FORMAT (8X,2A10,15X,I7,A10)
 420  FORMAT (49X,F11.2)
 430  FORMAT (49X,F11.6)
 440  FORMAT (8X,F20.12,6X,I3,1X,I2,1X,I2,I3,1X,I2,1X,F10.7)
      RETURN
      END
C*****************************************************************
      SUBROUTINE CRAMER(A,B,X)
C     
C     PURPOSE: SOLVE 3X3 SYSTEM A.X = B USING CRAMER'S RULE
C     
C     CODED BY: GERARD L.H. KRUIZINGA    06/15/1993
C     CENTER FOR SPACE RESEARCH    UNIVERSITY OF TEXAS AT AUSTIN
C     
C     INPUT:
C     
C     A         3X3 MATRIX
C     B         3X1 VECTOR
C     
C     OUTPUT:
C
C     X         3X1 SOLUTION  X = A-1.B
C
      IMPLICIT NONE

      DOUBLE PRECISION A(3,3),HELP(3,3),X(3),B(3)
      DOUBLE PRECISION DETA,DETX,DETY,DETZ

      INTEGER       I,J

C...  CHECK IF SOLUTION EXITS

      CALL DETERMINANT(A,DETA)

      IF (DETA .EQ. 0.0) THEN
        WRITE(*,*) ' SOLUTION DOES N0T EXIST BECAUSE '
        WRITE(*,*) ' DETERMINANT OF : '
        CALL WRITEMAT(A,3)
        WRITE(*,*) ' IS ZERO'
        RETURN
      ENDIF

      DO 20 I=1,3
        DO 10 J=1,3 
          HELP(I,J)=A(I,J)
          IF (J .EQ. 1) HELP(I,J)=B(I)
 10     CONTINUE
 20   CONTINUE

      CALL DETERMINANT(HELP,DETX)

      X(1) = DETX/DETA

      DO 40 I=1,3
        DO 30 J=1,3 
          HELP(I,J)=A(I,J)
          IF (J .EQ. 2) HELP(I,J)=B(I)
 30     CONTINUE
 40   CONTINUE

      CALL DETERMINANT(HELP,DETY)

      X(2) = DETY/DETA

      DO 60 I=1,3
        DO 50 J=1,3 
          HELP(I,J)=A(I,J)
          IF (J .EQ. 3) HELP(I,J)=B(I)
 50     CONTINUE
 60   CONTINUE

      CALL DETERMINANT(HELP,DETZ)

      X(3) = DETZ/DETA

      RETURN
      END
C*****************************************************************
      SUBROUTINE DETERMINANT(A,DET)
C     
C     PURPOSE: DETERMINE DETERMINANT OF A 3X3 MATRIX
C     
C     CODED BY: GERARD L.H. KRUIZINGA    06/15/1993
C     CENTER FOR SPACE RESEARCH    UNIVERSITY OF TEXAS AT AUSTIN
C     
C     INPUT:
C     
C     A         3X3 MATRIX
C     
C     OUTPUT:
C
C     DET       DETERMINANT OF MATRIX A
C
      IMPLICIT NONE

      DOUBLE PRECISION A(3,3),DET

      DET = A(1,1)*A(2,2)*A(3,3) +
     .      A(1,2)*A(2,3)*A(3,1) +
     .      A(1,3)*A(2,1)*A(3,2) -
     .      A(1,3)*A(2,2)*A(3,1) -
     .      A(1,1)*A(2,3)*A(3,2) -
     .      A(1,2)*A(2,1)*A(3,3)
 
      RETURN
      END
C*****************************************************************
      SUBROUTINE WRITEMAT (A,N)
C     
C     PURPOSE: WRITE OUT A 3X3 MATRIX A
C     
C     CODED BY: GERARD L.H. KRUIZINGA    06/15/1993
C     CENTER FOR SPACE RESEARCH    UNIVERSITY OF TEXAS AT AUSTIN
C
      IMPLICIT NONE
   
      DOUBLE PRECISION A(N,N)

      INTEGER       N,IROW,JCOL

      CHARACTER*1   NUMB(0:9)
      CHARACTER*80  FORMAT

      DATA NUMB /'0','1','2','3','4','5','6','7','8','9'/

      IF (N .GT. 9) THEN
        WRITE(*,*) ' SUBROUTINE WRITEMAT CAN ONLY PRINT ARRAYS'
        WRITE(*,*) ' WITH DIMENSION LESS OR EQUAL THAN 9'
        WRITE(*,*) ' CURRENT DIMENSION : ',N
        RETURN
      ENDIF

      FORMAT(1:1)='('
      FORMAT(2:2)=NUMB(N)
      FORMAT(3:7)='F9.3)'

      WRITE(6,*) ' '
      DO 10 IROW=1,N
        WRITE(6,FORMAT)(A(IROW,JCOL),JCOL=1,N)
 10   CONTINUE
      WRITE(6,*) ' '

      RETURN
      END
C*****************************************************************
      SUBROUTINE LEAPSEC(TREC,TORB,EPOCHREC,DY2,DY3,DY4,
     &                   A,B,C,FLAG,DTNEW)
C     
C     PURPOSE: DETERMINE IF INPUT TIME IS CLOSE TO A TIME
C              WHEN A LEAPSECOND WAS INSERTED 
C     
C     CODED BY: GERARD L.H. KRUIZINGA    06/16/1993
C     CENTER FOR SPACE RESEARCH     
C     
C     NOTE:    LEAP SECONDS UP TO DATE UNTIL JULY 1993
C     
      IMPLICIT NONE

      INTEGER          NLEAP

      PARAMETER        (NLEAP=19)

      DOUBLE PRECISION TREC,TORB,EPOCHREC,MJDLIM(NLEAP),TIME
      DOUBLE PRECISION EPSLEAP,TIME1,TIME2,TPREV,DTNEW
      DOUBLE PRECISION TRECMJD

      DOUBLE PRECISION A,B,C,DY2,DY3,DY4,DT1,DT2,DT3
      DOUBLE PRECISION X(3),Y(3),MAT(3,3)

      INTEGER          INDEX

      LOGICAL          FLAG,FIRST

      DIMENSION        TORB(4)

      DATA MJDLIM     /41317.0E+0, 41499.0E+0, 41683.0E+0,
     .                 42048.0E+0, 42413.0E+0, 42778.0E+0, 43144.0E+0,
     .                 43509.0E+0, 43874.0E+0, 44239.0E+0, 44786.0E+0,
     .                 45151.0E+0, 45516.0E+0, 46247.0E+0, 47161.0E+0,
     .                 47892.0E+0, 48257.0E+0, 48804.0E+0, 49169.0E+0
     .                /
      DATA EPSLEAP    /8640.0D0/
      DATA INDEX      /0/
      DATA TPREV      /0.0D0/
      DATA FIRST      /.TRUE./

C     DETERMINE THE LEAP-SECOND ARRAY INDEX WHICH IS AHEAD FOR
C     THE INPUT TIME WHEN LEAPSEC IS CALLED FOR THE FIRST TIME.  
C     THIS MEANS THAT THE FIRST TIME OF THE INTERPOLATION WINDOW 
C     MUST BE .GE. TIME OF LEAP + 10.0 SEC IN UTC

      IF (FIRST .OR. (TPREV - TREC) .GT. 0) THEN
         TIME1 = TORB(1)/86400.0D0 + EPOCHREC
 5       IF (INDEX .GE. NLEAP) GOTO 10
         INDEX = INDEX + 1
         IF (TIME1 .GE. MJDLIM(INDEX) + 1.0/EPSLEAP) GOTO 5
 10      CONTINUE
         FIRST = .FALSE.
      ENDIF

      FLAG  = .FALSE.
      TPREV = TREC
      TIME2 = TORB(4)/86400.0D0 + EPOCHREC

C     CHECK IF LAST TIME OF INTERPOLATION WINDOW IS .LE. THAN
C     TIME OF LEAP-SECOND

      IF (TIME2 .LE. MJDLIM(INDEX)) RETURN

C     CHECK IF FIRST TIME OF INTERPOLATION WINDOW IS .GE. THAN
C     TIME OF LEAP-SECOND + 10.0  IN UTC
   
      TIME1 = TORB(1)/86400.0D0 + EPOCHREC
      IF (TIME1 .GE. MJDLIM(INDEX) + 1.0/EPSLEAP) THEN
        INDEX = INDEX + 1
        IF (INDEX .GT. NLEAP) INDEX = NLEAP
        RETURN
      ENDIF

C     INTERPOLATION AT THIS POINT CONTAINS A 11 SECOND INTERVAL 
C     NOTIFY USER OF OCCURENCE OF THE LEAP SECOND

      WRITE(*,*) ' '
      WRITE(*,*) ' LEAP SECOND ENCOUNTERED !!!!'
      WRITE(*,*) ' TIME OF LEAP-SECOND (MJD)  : ',MJDLIM(INDEX)
      WRITE(*,*) ' FIRST TIME OF INTERVAL(MJD): ',TIME1
      WRITE(*,*) ' LAST TIME OF INTERVAL(MJD) : ',TIME2
      WRITE(*,*) ' PLEASE VERIFY RESULTS !!!!!'
      WRITE(*,*) ' '
 
C     DETERMINE IN WHICH INTERVAL THE 11 SECONDS OCCURS

      TRECMJD = TREC/86400.0D0 + EPOCHREC

C     INTERVAL BETWEEN TORB(1),TORB(2)
      TIME = TORB(2)/86400.0D0 + EPOCHREC
      IF ( MJDLIM(INDEX) .LE. TIME) THEN
        DT1 = 11.0
        DT2 = 21.0
        DT3 = 31.0
        GOTO 20
      ENDIF

C     INTERVAL BETWEEN TORB(2),TORB(3)
      TIME = TORB(3)/86400.0D0 + EPOCHREC
      IF ( MJDLIM(INDEX) .LE. TIME) THEN
        DT1 = 10.0
        DT2 = 21.0
        DT3 = 31.0
        GOTO 20
      ENDIF

C     INTERVAL BETWEEN TORB(3),TORB(4)
      DT1 = 10.0
      DT2 = 20.0
      DT3 = 31.0
      GOTO 20

 20   CONTINUE

      IF (TRECMJD .GT. MJDLIM(INDEX)) THEN
        DTNEW = TREC - TORB(1) + 1.0
      ELSE
        DTNEW = TREC - TORB(1)
      ENDIF  

      MAT(1,1) = DT1*DT1*DT1
      MAT(2,1) = DT2*DT2*DT2
      MAT(3,1) = DT3*DT3*DT3

      MAT(1,2) = DT1*DT1
      MAT(2,2) = DT2*DT2
      MAT(3,2) = DT3*DT3

      MAT(1,3) = DT1
      MAT(2,3) = DT2
      MAT(3,3) = DT3

      Y(1)=DY2
      Y(2)=DY3
      Y(3)=DY4

      CALL CRAMER(MAT,Y,X)

      A = X(1)
      B = X(2)
      C = X(3)

      FLAG = .TRUE.
      RETURN
      END
