       SUBROUTINE FILMAIN(TRIPS,HEAD,FTRIPS,COFINM,IOS)
*
*************************************************************************
*	
*      TITLE:  FILTER MAIN
*
*      FILE NAME:  FILMAIN.FOR
*
*      PURPOSE:  MAIN PROGRAM THAT RECEIVES THE 1.92 SECOND AVERAGED
*                       DATA TO BE FILTERED.
*
*     HISTORY:
*
*        AUTHOR               DATE            CHANGE
*        ------               ----            ------
*        S.J.KEMPLER          12/1/88         VERSION 1.0; ORIGINAL CODE
*                             2/16/89         VERSION 2.0; NEW INTERPOLATION
*                                             METHOD; PARAMETER FOR DEGREE
*                                             OF LEAST SQUARES FIT ADDED;
*                                             VERSION NUMBER PRINTED OUT.
*                             3/1/89          VERSION 2.1: FIX TO CHECK
*                                             BEGIN OF BLOCK OF DATA WHEN
*                                             BLOCK ENDS WITH FILL-IN FINDFILL
*                             3/14/89         V.2.2: FIX TO CHECK FOR SMALL
*                                             GAPS AFTER CONTINUOUS DATA
*                                             BREAK; CHECK END-OF-DATA
*                                             FLAG BEFORE ENTERING FINDGAP
*                             4/10/89         V.2.3: PRINTOUT NUMBER OF
*                                             COEFFICIENTS.
*
*      CALLING SEQUENCE:  SUBROUTINE FILMAIN(TRIPS,HEAD,FTRIPS,COFINM,IOS)
*
*      ARGUMENT LIST:
*
*       NAME               TYPE   USE    PURPOSE
*       ----               ----   ---    -------
*
*	TRIPS(3,0:19999)   R*4      I    DATA TO BE FILTERED.
*
*	HEAD(64,0:799)     I*2      I    HEADER RECORD CORRESPONDING TO 25
*                                        DATA TRIPLETS IN TRIPS AND FTRIPS.
*
*	FTRIPS(3,0:19999)  R*4      O    FILTERED DATA.
*
*	COFINM             C*8      O    CONTAINS ALPHANUMERIC NAME OF
*                                        FILTER COEFFICIENTS FILE USED.
*       IOS                I*4      O    FILE OPEN STATUS. IOS NOT EQUAL
*                                        TO 0, WILL GENERATE REPORT
*                                        IMMEDIATELY AND TERMINATE PROGRAM.
*
*      MODULES REFERENCED:  REDFIL, FINDGAPS, GENREP
*
*      NOTES:  THIS ROUTINE RECIEVES THE BUFFER OF 20000 TRIPLETS OF DATA
*                     AND PASSES THEM TO SUBSEQUENT PROGRAMS TO BE
*                     INTERPOLATED AND FILTERED.  TH E PROCESSED TRIPLETS ARE
*                     PASSED BACK IN A DUPLICATE BUFFER CALLED FTRIPS.  THE
*                     FIRST TIME CALLED, 10000 PLUS THE NUMBER OF FILTER
*                     COEFFICIENTS OF DATA ARE PROCESSED AND PASSED BACK.
*                     ALL SUBSEQUENT CALLS WILL PROCESS THE NEXT 10000
*                     TRIPLETS.
*
*
*      PDL:
*
*       INITIALIZE STATS TO 0
*       IF THIS IS THE FIRST ACCESS OF FILMAIN
*       THEN
*          CALL REDFIL TO READ PARAMETER FILES
*       ENDIF
*       IF IOS EQUALS 0
*       THEN
*          CALL FINDGAPS TO INTERPOLATE AND FILTER DATA
*       ENDIF
*       IF END OF DATA IS REACHED OR IOS NOT EQUAL TO 0
*       THEN
*          CALL GENREP TO GENERATE RUN STATISTICS REPORT
*       ENDIF
*       RETURN
*
*******************************************************************
C
      REAL*4 TRIPS(3,0:19999), FTRIPS(3,0:19999), COEF(500)
      INTEGER*2 HEAD(64,0:799)
      CHARACTER*8 COFINM
      INTEGER*4 STATS(10),NUMCOF,DATGAP,INTPTS,FITDEG,EODFLG/0/
C
      DATA IFIRST/0/
C
      IOS=0
C
C      PRINT *,' IN FILTMAIN:'
C      PRINT *,(JJ,FTRIPS(1,JJ),JJ=0,19999,20)
C
C     READ INPUT FILES FIRST TIME ONLY
C
      IF (IFIRST.EQ.0) THEN
          IFIRST=1
          CALL REDFIL (HEAD,COEF,NUMCOF,DATGAP,INTPTS,COFINM,FITDEG,
     C                 IOS)
C
          DO 5 I=1,10
              STATS(I)=0
  5       CONTINUE
C
      PRINT *,' ***************************************************'
      PRINT *,' HERE COMES THE VOYAGER DIGITAL FILTERING PROGRAM!!'
C
      ENDIF
C
C     PRINT *,' EODFLG,IOS: ', EODFLG,IOS
      IF(EODFLG.NE.0) GO TO 22
C
C     IF(IOS.NE.0) GO TO 20
C
c      WRITE(98,*) ' BEFORE FINDGAPS'
c      WRITE(98,*) ' I,TRIPS(1:3,I),FTRIPS(1:3,I):'
c      DO I = 1,10000
c       K = (I-1)/25 + 1
c       WRITE(98,'(3(1X,I5),2(3(1X,F7.3),3X))') 
c     P           I,HEAD(25,K),HEAD(26,K),
c     P           (TRIPS(J,I),J=1,3),(FTRIPS(J,I),J=1,3)
c      END DO
      CALL FINDGAPS(TRIPS,HEAD,FTRIPS,COEF,NUMCOF,DATGAP,INTPTS,
     C              STATS,EODFLG,FITDEG)
C
C     WRITE(99,*) ' AFTER FINDGAPS:'
C     WRITE(99,*) ' I,TRIPS(1:3,I),FTRIPS(1:3,I):'
C     DO I = 1,10000
C      K = (I-1)/25 + 1
C      WRITE(99,'(3(1X,I5),2(3(1X,F7.3),3X))') 
C    P           I,HEAD(25,K),HEAD(26,K),
C    P           (TRIPS(J,I),J=1,3),(FTRIPS(J,I),J=1,3)
C     END DO
C      PRINT *,' TRIPS(Y,10001-10100),FTRIPS(Y,10001-10100):',
C     P         (I,TRIPS(2,I),FTRIPS(2,I),I=10001,10100)
C      PRINT *,' STATS:',STATS
C
C     OUTPUT PROGRAM STATISTICS LAST TIME ONLY
C
C 20   IF(EODFLG.EQ.1.OR.IOS.NE.0) CALL GENREP(STATS,COFINM,NUMCOF)
 20   IF(EODFLG.EQ.1) CALL GENREP(STATS,COFINM,NUMCOF)
C
 22   RETURN
      END
C
C
C
C
C
C
      SUBROUTINE REDFIL(HEAD,COEF,NUMCOF,DATGAP,INTPTS,COFINM,
     S                  FITDEG,IOS)
*************************************************************************
* 	
*      TITLE:  READ DATA FILES
*
*      FILE NAME:  REDFIL.FOR
*
*      PURPOSE:  READS INPUT PARAMETER FILES AND PASSES THEM ON.
*
*      HISTORY:
*
*        AUTHOR               DATE            CHANGE
*        ------               ----            ------
*        S.J.KEMPLER          12/1/88         ORIGINAL CODE
*                             2/16/89         ADD DEGREE OF FIT PARAMETER
*
*      CALLING SEQUENCE:  SUBROUTINE REDFIL(HEAD,COEF,NUMCOF,DATGAP,
*                             INTPTS,COFINM,FITDEG,IOS)
*
*      ARGUMENT LIST:
*
*       NAME               TYPE   USE    PURPOSE
*       ----               ----   ---    -------
*	HEAD(64,0:799)     I*2     I     HEADER RECORD CORRESPONDING TO 25
*                                        DATA TRIPLETS IN TRIPS AND FTRIPS.
*       COEF(500)          R*4     O     THE ENTIRE ARRAY OF FILTER COEF TO BE
*                                        USED (2N-1, IF SYMMETRICAL; N IF NOT
*                                        SYMMETRICAL, WHERE N IS THE NUMBER
*                                        OF COEFF RECIEVED FROM THE FILTER
*                                        COEF FILE.  ONLY HALF OF THE COEF
*                                        ARE PROVIDED IF THE FILTER IS
*                                        SYMMETRICAL.)
*        NUMCOF            I*4     O     TOTAL NUMBER OF FILTER COEF (SEE
*                                        ABOVE).
*        DATGAP            I*4     O     THE NUMBER OF CONSECUTIVE FILL DATA
*                                        VALUES (999.0) PERMITTED BEFORE THE
*                                        "GAP" IS TOO LARGE TO INTERPOLATE.
*                                        (24 MAX)
*        INTPTS            I*4     O     THE NUMBER OF ADJACENT VALID DATA
*                                        POINTS TO USE ON EITHER SIDE OF A
*                                        PERIOD
*        COFINM            C*8     O     CONTAINS ALPHANUMERIC NAME OF
*                                        FILTER COEFFICIENTS FILE USED.
*        FITDEG            I*4     O     DEGREE OF LEAST SQUARES FIT
*        IOS               I*4     O     FILE OPEN STATUS. IOS NOT EQUAL
*                                        TO 0, WILL GENERATE REPORT
*                                        IMMEDIATELY AND TERMINATE PROGRAM.
*
*      MODULES REFERENCED: NONE
*
*      DATA FILES ACCESSED:
*
*        NAME            PURPOSE
*        ----            -------
*        COEFNA.DAT      CONTAINS A LIST OF FILTER COEFFICIENT FILE NAMES
*                          AND THE END DATE ON WHICH THE COEFFICIENTS
*                          ARE VALID;  ONE PER RECORD IN INCREASING
*                          SEQUENTIAL ORDER BY DATE.  EACH RECORD CONTAINS:
*                 IYEAR, IDAY, COFINM
*                          THESE PARAMETERS ARE READ INTO THE PROGRAM WITH
*                          THE FOLLOWING FORMAT:
*                 I2, 1X, I3, A8
*
*        THE FILTER COEFFICIENT FILE NAME PROVIDED BY COEFNA.DAT.  THIS FILE
*                                 CONTAINS:
*          RECORD 1:  ISYMFL - (0=ASSYMMERICAL FILTER COEF; 1=SYMMETRICAL)
*                     NUMCOF - NUMBER OF FILTER COEFFICIENTS TO FOLLOW.
*                              RECORD 1 FORMAT: I1,1X,I3
*          RECORDS 2 TO NUMCOF+1 - ONE PER LINE, THE FILTER COEF PROVIDED. READ
*                          IN WITH AN F15.10 FORMAT (READS ONE FLOATING POINT).
*
*        PROPAR.DAT      CONTAINS  VARIOUS PROGRAM PARAMETERS, ONE PER RECORD:
*                DATGAP - THE NUMBER OF CONSECUTIVE FILL DATA VALUES
*                         (999.0) PERMITTED BEFORE THE "GAP" IS TOO
*                         LARGE TO INTERPOLATE. (24 MAX.).  READ WITH I2 FORMAT
*                INTPTS - THE NUMBER OF ADJACENT VALID DATA POINTS
*                         TO USE ON EITHER SIDE OF A PERIOD OF DATA TO
*                         BE LEAST SQUARES INTERPOLATED.  READ WITH I2 FORMAT.
*                FITDEG - DEGREE OF LEAST SQUARES FIT (5 MAX). I2 FORMAT.
*
*
*      NOTES:  COEFNA.DAT IS READ FIRST TO RECIEVE THE NAME OF THE DESIRED
*           FILTER COEFFICIENT FILE.  THEN THE FILTER COEFFICIENT FILE IS READ.
*           THE NUMBER OF COEFFICIENTS RECIEVED REPRESENTS ONE HALF
*           OF A SYMMETRICAL FILTER IF SYMMETRY FLAG IS FLAG IS SET OR
*           ALL THE COEF OF A FILTER IF FLAG IS NOT SET.  IF THE FORMER IS TRUE,
*           THE ENTIRE FILTER ARRAY IS PROVIDED.
*
*      PDL:
*         INITIALIZE IYEAR AND IDAY TO 0
*         DO UNTIL DATE IN COEFNA.DAT IS GREATER THAN FIRST DATE IN HEAD
*             READ COEFNA.DAT FILE
*         ENDDO
*         READ FILTER COEF FILE
*         IF ISYMFL EQUALS 1
*         THEN
*             DO I EQUAL 1 TO (NUMCOF-1)
*                 SET COEF(NUMCOF+I) TO COEF(NUMCOF-I)
*             ENDDO
*             SET NUMCOF EQUAL TO (2 * NUMCOF - 1)
*         ENDIF
*         READ PROPAR.DAT FILE
*         RETURN
*******************************************************************
C
      REAL*4 COEF(500)
      INTEGER*4 NUMCOF,DATGAP,INTPTS,ISYMFL,FITDEG
      INTEGER*2 STYEAR,STDATE
      INTEGER*2 HEAD(64,0:799)
      CHARACTER*8 COFINM
C
C     RETRIEVE DATE AND TIME FIRST OF FIRST RECORD AND SEARCH FOR COEFICIENTS
C     FILE NAME THAT CORRESPONDS TO THAT TIME PERIOD
C
      STYEAR=HEAD(7,0)
      STDATE=HEAD(8,0)
      ISTY4  =STYEAR
      ISTD4  =STDATE
      INDATE = (ISTY4*1000)+ISTD4
C
C   READ FILTER COEFFICIENTS FILE
C
 101  FORMAT(I1,1X,I3)
      READ(9,101) ISYMFL,NUMCOF
      PRINT *, 'ISYMFL,NUMCOF',ISYMFL,NUMCOF
C
 102  FORMAT(F15.10)
      DO 10 I=1,NUMCOF
          READ(9,102) COEF(I)
C         PRINT *, 'I,COEF(I)',I,COEF(I)
  10  CONTINUE
C
C     IF DATA IS TO BE MADE SYMMETRICAL...
C
      IF(ISYMFL.EQ.1) THEN
          DO 15 I=1,NUMCOF-1
              COEF(NUMCOF+I) = COEF(NUMCOF-I)
  15      CONTINUE
          NUMCOF = (2*NUMCOF) -1
      ENDIF
C
C   READ PROGRAM PARAMETERS FILE
C
      READ(10,103) DATGAP
      READ(10,103) INTPTS
      READ(10,103) FITDEG
C     PRINT *, 'DATGAP,INTPTS,FITDEG',DATGAP,INTPTS,FITDEG
 103  FORMAT(I2)
C     PRINT *,DATGAP,INTPTS,ISYMFL,NUMCOF,(COEF(I), I=1,NUMCOF+2)
C
      RETURN
      END
C
C
C
C
C
      SUBROUTINE FINDGAPS(TRIPS,HEAD,FTRIPS,COEF,NUMCOF,
     S                 DATGAP,INTPTS,ISTATS,EODFLG,FITDEG)
*************************************************************************
*	
*      TITLE:  FIND DATA GAPS
*
*      FILE NAME:  FINDGAPS.FOR
*
*      PURPOSE:  SEARCHES THE BLOCK OF DATA FOR MISSING DAT - GAPS OF 25
*                          TRIPLETS IN LENGTH.
*
*      HISTORY:
*
*        AUTHOR               DATE            CHANGE
*        ------               ----            ------
*        S.J.KEMPLER          12/1/88         ORIGINAL CODE
*                             2/16/89         ADDED FITDEG TO ARGUMENT LIST
*        T.P.MCCLANNAHAN      3/08/94         REMOVE I*2/I*4 EQUIVALENCE AND
*        S.B.KRAMER                           REPLACE WITH DIRECT ASSIGNMENT
*
*      CALLING SEQUENCE:  SUBROUTINE FINDGAPS(TRIPS,HEAD,FTRIPS,COEF,NUMCOF
*                           DATGAP,INTPTS,ISTATS,EODFLG,FITDEG)
*
*      ARGUMENT LIST:
*
*       NAME               TYPE   USE    PURPOSE
*       ----               ----   ---    -------
*
*	TRIPS(3,0:19999)   R*4     I     DATA TO BE FILTERED.
*
*	HEAD(64,0:799)     I*2     I     HEADER RECORD CORRESPONDING TO 25
*                                          DATA TRIPLETS IN TRIPS AND FTRIPS.
*
*	FTRIPS(3,0:19999)  R*4     O     FILTERED DATA.
*       COEF(500)          R*4     O     THE ENTIRE ARRAY OF FILTER COEF TO BE
*                                        USED. (2N-1, IF SYMMETRICAL; N IF NOT
*                                        SYMMETRICAL, WHERE N IS THE NUMBER
*                                        OF COEFF RECIEVED FROM THE FILTER COEF
*                                        FILE.  ONLY HALF OF THE COEF ARE PRO-
*                                        VIDED IF THE FILTER IS SYMMETRICAL.)
*        NUMCOF            I*4     I     TOTAL NUMBER OF FILTER COEF (SEE
*                                        ABOVE).
*        DATGAP            I*4     I     THE NUMBER OF CONSECUTIVE FILL DATA
*                                        VALUES (999.0) PERMITTED BEFORE THE
*                                        "GAP" IS TOO LARGE TO INTERPOLATE.
*                                        (24 MAX)
*        INTPTS            I*4     I     THE NUMBER OF ADJACENT VALID DATA
*                                        POINTS TO USE ON EITHER SIDE OF A
*                                        PERIOD OF DATA TO BE SPLINE
*                                        INTERPOLATED.
*        ISTATS(10)        I*4     O     ISTATS(1)=NUMBER OF INTERPOLATED
*                                        TRIPLETS ADDED; ISTATS(2)=NUMBER OF
*                                        INTERPOLATED PERIODS; ISTATS(3)=N0.OF
*                                        FILL TRIPLETS ADDED; ISTATS(4)= NUMBER
*                                        OF GAPS AND LARGE MISSING DATA PERIODS.
*                                        ISTATS(5)=TOTAL NUMBER OF TRIPS.
*        EODFLG            I*4     O     FLAG=1 INDICATES END OF DATA
*        FITDEG            I*4     I     DEGREE OF LEAST SQUARES FIT
*
*      MODULES REFERENCED:  FINDFILL
*
*      NOTES:   THE TWO FDS COUNTERS ARE PULLED OFF EACH HEADER AND COMPARED
*                    TO PROGRAM COUNTER TO SEARCH FOR GAPS.  IF FOUND, PROGRAM
*                    CONTROL GOES TO SUBROUTINE FINDGAPS FOR FURTHER PROCESSING.
*
*      PDL:
*         IF FIRST TIME CALLED
*         THEN
*             SET ITERATE TO 400 + (NUMBER OF FILTER COEFFICIENTS/25+1) + 1
*             INITILAIZE COUNT TO 0
*             INITIALIZE NUMTRP TO 0
*             INITIALIZE EODFLG TO 0
*         ELSE
*             INITIALIZE NUMTRP TO 0
*             SET ITERATE TO 400
*         ENDIF
*         INITIALIZE REPEAT TO 1
*         SET ICNT1 TO THE FASTER INCREMENTING FDS COUNTER IN THE NEXT
*           HEADER RECORDTO BE ACCESSED
*         SET ICNT2 TO THE SLOWER INCREMENTING FDS COUNTER IN THE NEXT
*           HEADER RECORDTO BE ACCESSED
*         DO WHILE REPEAT IS LESS THAN 50000
*             INCREMENT REPEAT BY 1
*             INCREMENT NUMTRP BY 25
*             INCREMENT ICNT1 BY 1
*             IF ICNT1 IS GREATER THAN 59
*             THEN
*                 SET ICNT1 TO 0
*                 INCREMENT ICNT2 BY 1
*             ENDIF
*             IF (FDS COUNTERS EQUAL 0)
*             THEN
*                CALL FINDFILL TO SEARCH FOR FILL DATA
*                SET STAT(5) EQUAL TOP NUMTRP
*                SET EODFLG TO 1
*                RETURN
*             ENDIF
*             IF (ICNT1 OR INCT2 DO NOT EQUAL THE FDS COUNTERS IN THE NEXT
*               HEADER RECORD) OR (REPEAT EQUALS ITERATE)
*             THEN
*                 IF (ICNT1 AND INCT2 ARE EQUAL TO THE FDS COUNTERS IN THE NEXT
*                   HEADER RECORD)
*                 THEN
*                     IF(LAST TRIPLET TO BE SENT IS FILL)
*                     THEN
*                        INCREMENT ITERATE BY 1
*                        GO TO BEGINNING OF LOOP
*                     ENDIF
*                     SET CONTINUATION FLAG TO 1
*                 ELSE
*                     SET CONTINUATION FLAG TO 0
*                 ENDIF
*                 CALL FINDFILL TO SEARCH FOR FILL DATA
*                 SET COUNT TO COUNT + NUMTRP
*                 RESET NUMTRP TO 0
*                 SET ICNT1 TO THE FASTER INCREMENTING FDS COUNTER IN THE LAST
*                   HEADER RECORD ACCESSED
*                 SET ICNT2 TO THE SLOWER INCREMENTING FDS COUNTER IN THE LAST
*                   HEADER RECORD ACCESSED
*                 IF (FDS COUNTERS EQUAL 0) OR (REPEAT EQUALS ITERATE)
*                 THEN
*                     RETURN
*                 ENDIF
*             ENDIF
*         ENDDO
*         RETURN
*******************************************************************
C
      REAL*4 TRIPS(3,0:19999),FTRIPS(3,0:19999),COEF(500)
      INTEGER*2 HEAD(64,0:799)
      INTEGER*4 NUMCOF,DATGAP,INTPTS,ISTATS(10),EODFLG,FITDEG
      INTEGER*4 ITERAT,COUNT,NUMTRP,REPEAT,ICNT2S/0/,SLONX4/0/
      INTEGER*2 ICNT1F,SLOTMP(2),SLONXT(2),STDATE,STHOUR,STMIN
C
C REMOVED AND REPLACED WITH DIRECT ASSIGNMENTS  03/08/94 TPM
C
C     EQUIVALENCE(ICNT2S,SLOTMP(1))
C     EQUIVALENCE(SLONX4,SLONXT(1))
C
      DATA IFIRST/0/
      DATA FILL/999.0/
C
      PRINT *,' IN FINDGAPS:'
      IF (IFIRST.EQ.0) THEN
          IFIRST=1
          EODFLG = 0
          COUNT=0
          NUMTRP= 0
          ITERAT= 400 + (NUMCOF/25) + 2
          REPEAT=0
          ICNT1F=HEAD(26,0)
          SLOTMP(2)=HEAD(25,0)
          ICNT2S = SLOTMP(2)  ! ADDED 03/08/94 TPM
C
          IEXTRA = 0
C
      ELSE
          NUMTRP= 0
          ITERAT= ITERAT + 400 - IEXTRA
          IEXTRA = 0
      ENDIF
C
C     REPEAT FOR THE NUMBER OF RECORDS TO BE SEARCHED IN THIS CALL
C
      DO  5 I=1,50000
C
          REPEAT = REPEAT + 1
          SLONXT(2) = HEAD(25,MOD(REPEAT,800))
          SLONX4 = SLONXT(2)  ! ADDED 03/08/94 TPM
C
          NUMTRP = NUMTRP + 25
          ICNT1F = ICNT1F + 1
          IF (ICNT1F.GT.59) THEN
              ICNT1F = 0
              ICNT2S = ICNT2S + 1
          ENDIF
C
C     IF BOTH RECORD COUNTERS EQUAL ZERO, END OF DATA REACHED. GO SEARCH
C     FOR MISSING DATA, THEN RETURN.
C
          IF (HEAD(26,MOD(REPEAT,800)).EQ.0.AND.
     I        SLONX4.EQ.0) THEN
              ICONT = 0
              CALL FINDFILL(TRIPS,FTRIPS,COEF,NUMCOF,COUNT,NUMTRP,
     C                       DATGAP,INTPTS,ISTATS,ICONT,FITDEG)
              EODFLG = 1
              ISTATS(5) = COUNT + NUMTRP
              RETURN
          ENDIF
C
C     IF NEITHER OF THE COUNTERS SKIP AND WE HAVE NOT REACHED THE LAST
C     ITERATION OF THIS CALL, SKIP DOWN AND LOOK AT THE NEXT DATA RECORD
C
          IF (ICNT1F.EQ.HEAD(26,MOD(REPEAT,800)).
     I   AND .ICNT2S.EQ.SLONX4.
     I   AND .REPEAT.LT.ITERAT) GO TO 40
C
C     GET DATE AT BEGINNING OF CURRENT BLOCK OF DATA TO BE FILTERED
C
          STDATE = HEAD(8, MOD(REPEAT,800))
          STHOUR = HEAD(9, MOD(REPEAT,800))
          STMIN  = HEAD(10,MOD(REPEAT,800))
C
          PRINT *,' ** DATE= ',STDATE,' HOUR= ',STHOUR,' MIN= ',STMIN
C
      PRINT *,' ** REPEAT,FAST COUNTER,MY FAST COUNTER,SLOW,MY SLOW:'
      PRINT *,REPEAT,HEAD(26,MOD(REPEAT,800)),ICNT1F,SLONX4,ICNT2S
C
C     IF THE NEXT STATEMENT TRUE, MUST BE AT LAST ITERATION.
C
              IF (ICNT1F.EQ.HEAD(26,MOD(REPEAT,800)).AND.ICNT2S.EQ.
     I                      SLONX4 ) THEN
C
               LAST = COUNT +NUMTRP -1
C
      PRINT *,' ITERAT,REPEAT,TRIPS(2,LAST),LAST,COUNT,
     P          NUMTRP,I,ICONT:'
      PRINT *,ITERAT,REPEAT,TRIPS(2,MOD(LAST,20000)),LAST,
     P        COUNT,NUMTRP,I,ICONT
C
C     JUST IN CASE THE LAST TRIPLET IS FILL, GO GET ANOTHER RECORD...
C     SUBROUTINE FINDFILL NEEDS GOOD DATA AFTER FINDING FILL
C
               IF (TRIPS(1,MOD(LAST,20000)).EQ.FILL. OR.
     I             TRIPS(2,MOD(LAST,20000)).EQ.FILL. OR.
     I             TRIPS(3,MOD(LAST,20000)).EQ.FILL) THEN
                  ITERAT = ITERAT + 1
                  IEXTRA = IEXTRA + 1
C
      PRINT *,' ITERAT,REPEAT,TRIPS(2,LAST),LAST,COUNT,
     P          NUMTRP,I,ICONT:'
      PRINT *,ITERAT,REPEAT,TRIPS(2,MOD(LAST,20000)),LAST,
     P        COUNT,NUMTRP,I,ICONT
                  GO TO 40
C
                ENDIF
                  ICONT = 1
              ELSE
                  ICONT = 0
              ENDIF
      PRINT *,' COUNT,NUMTRP,HEAD(26),HEAD(25),REPEAT,ITERAT:'
      PRINT *, COUNT,NUMTRP,HEAD(26,MOD(REPEAT,800))
      PRINT *, SLONX4,REPEAT,ITERAT
C
C      CALL FINDFILL TO SEARCH THE DATA ALREADY EXAMINED FOR GAPS, FOR FILL
C
              CALL FINDFILL(TRIPS,FTRIPS,COEF,NUMCOF,COUNT,NUMTRP,
     C                       DATGAP,INTPTS,ISTATS,ICONT,FITDEG)
C
C     UPDATE POINTERS...IF LAST ITERATION RETURN
C
              COUNT = COUNT + NUMTRP
              NUMTRP = 0
              ICNT1F=HEAD(26,MOD(REPEAT,800))
              ICNT2S=SLONX4
              IF (REPEAT.EQ.ITERAT) RETURN
C
  40      CONTINUE
  5   CONTINUE
      RETURN
      END
C
C
C
C
C
      SUBROUTINE FINDFILL (TRIPS, FTRIPS, COEF, NUMCOF,
     S           ISTART, INUM, DATGAP, INTPTS, ISTATS, ICONT, FITDEG)
C
*************************************************************************
*	
*      TITLE:  FIND FILL DATA
*
*      FILE NAME:  FINDFILL.FOR
*
*      PURPOSE:  SEARCHES SECTIONS OF DATA FOR FILL DATA.
*
*      HISTORY:
*
*        AUTHOR               DATE            CHANGE
*        ------               ----            ------
*        S.J.KEMPLER          12/1/88         ORIGINAL CODE
*                             3/1/89          VERSION 2.1: FIX TO CHECK
*                                             BEGIN OF BLOCK OF DATA WHEN
*                                             BLOCK ENDS WITH FILL
*                             3/14/88         V.2.2: FIX TO CHECK FOR SMALL
*                                             GAPS AFTER CONTINUOUS DATA
*                                             BREAK; CHECK END-OF-DATA
*                                             FLAG BEFORE ENTERING FINDGAP
*
*      CALLING SEQUENCE:  SUBROUTINE FINDFILL (TRIPS, FTRIPS, COEF, NUMCOF,
*                     ISTART, INUM, DATGAP, INTPTS, ISTATS, ICONT, FITDEG)
*
*      ARGUMENT LIST:
*
*       NAME               TYPE   USE    PURPOSE
*       ----               ----   ---    -------
*
*	TRIPS(3,0:19999)   R*4     I     DATA TO BE FILTERED.
*	FTRIPS(3,0:19999)  R*4     O     FILTERED DATA.
*       COEF(500)          R*4     O     THE ENTIRE ARRAY OF FILTER COEF TO BE
*                                        USED. (2N-1, IF SYMMETRICAL; N IF NOT
*                                        SYMMETRICAL, WHERE N IS THE NUMBER
*                                        OF COEFF RECIEVED FROM THE FILTER COEF
*                                        FILE.  ONLY HALF OF THE COEF ARE PRO-
*                                        VIDED IF THE FILTER IS SYMMETRICAL.)
*       NUMCOF             I*4     I     TOTAL NUMBER OF FILTER COEF (SEE
*                                        ABOVE).
*       ISTART             I*4     I     THE FIRST TRIPLET IN TRIPS TO START
*                                        CHECKING FOR FILL DATA.
*       INUM               I*4     I     NUMBER OF TRIPLETS IN THIS GROUP TO BE
*                                        SEARCHED FOR FILL DATA.
*       DATGAP             I*4     I     THE NUMBER OF CONSECUTIVE FILL DATA
*                                        VALUES (999.0) PERMITTED BEFORE THE
*                                        "GAP" IS TOO LARGE TO INTERPOLATE.
*                                        (24 MAX)
*       INTPTS             I*4     I     THE NUMBER OF ADJACENT VALID DATA
*                                        POINTS TO USE ON EITHER SIDE OF A
*                                        PERIOD OF DATA TO BE SPLINE
*                                        INTERPOLATED.
*       ISTATS(10)         I*4     O     ISTATS(1)=NUMBER OF INTERPOLATED
*                                        TRIPLETS ADDED; ISTATS(2)=NUMBER OF
*                                        INTERPOLATED PERIODS; ISTATS(3)=N0.OF
*                                        FILL TRIPLETS ADDED; ISTATS(4)= NUMBER
*                                        OF GAPS AND LARGE MISSING DATA
*                                        PERIODS.
*                                        ISTATS(5)=TOTAL NUMBER OF TRIPS.
*       ICONT              I*4     I     CONTINUATION FLAG: 0=NEXT GROUP
*                                        OF TRIPLETS FOLLOWS DATA GAP; 1=
*                                        NEXT GROUP OF TRIPLETS IS A
*                                        CONTINUATION OF GOOD DATA FROM
*                                        LAST GROUP.
*        FITDEG            I*4     I     DEGREE OF LEAST SQUARES FIT
*
*      MODULES REFERENCED:  INTERP, FILTDAT
*
*      NOTES:  FINDFILL RECEIVES THE START TRIPLET AND THE NUMBER OF TRIPLETS TO
*           PROCESS IN THIS PARTICULAR CALL.  FILL DATA IS SEARCHED FOR.  ONCE
*           FOUND,  THE NUMBER OF CONSECUTIVE FILL TRIPLETS ARE RECORDED.
*           IF MORE THAN 24 TRIPLETS ARE FILL, CONTROL IS SENT TO FILTER ALL
*           DATA UP TO THIS TRIPLET.  IF NOT, THE FILL DATA IS INTERPOLATED FOR.
*
*      PDL:
*
*         SET JCONT TO 0
*         SET FLTBEG TO ISTART
*         SET FLTCNT TO 0
*         SET COUNTR TO ISTART
*         SET INTCNT TO 0
*         DO WHILE LESS THAN OR EQUAL TO INUM
*             INITIALIZE FTRIPS(ITERATION) TO FILL
*             INCREMENT FLTCNT BY 1
*             IF TRIPS(COUNTR) EQUALS FILL OR ITERATION INUM IS REACHED
*             THEN
*                 IF ITERATION INUM IS REACHED AND ICONT EQUALS 1
*                 THEN
*                     SET JCONT TO 1
*                 ELSE
*                     SET JCONT TO 0
*                 ENDIF
*                 IF INTCNT EQUALS 0
*                 THEN
*                     SET INTBEG TO COUNTR
*                 ENDIF
*                 INCREMENT INTCNT BY 1
*             ENDIF
*             IF JCONT EQUALS 0 AND INUM IS REACHED AND TRIPS EQUALS FILL
*             THEN
*                 SET FLTCNT TO FLTCNT - INTCNT
*                 CALL FILTDAT TO FILTER TRIPLETS SINCE LAST CALL UP TO
*                   PRESENT FILL GAP
*                 IF JCONT EQUALS 0
*                 THEN
*                     INCREMENT ISTATS(4) BY 1
*                 ENDIF
*                 SET INTCNT TO 0
*                 SET FLTBEG TO COUNTR
*                 SET FLTCNT TO 1
*             ENDIF
*             IF INTCNT IS GREATER THAN 0 AND TRIPS(COUNTR) IS NOT FILL
*             THEN
*                 DECREMENT INTCNT BY 1
*                 IF INTCNT IS GREATER THAN DATGAP OR ITERATION INUM IS REACHED
*                 THEN
*                     IF ITERATION INUM IS NOT REACHED
*                     THEN
*                         DECREMENT FLTCNT BY 1
*                     ENDIF
*                     SET FLTCNT TO FLTCNT - INTCNT
*                     IF FIRST FILTERING FOR BLOCK AND (LAST BLOCK
*                       IS SMALLER THAN (NUMCOF-1)/2)
*                     THEN
*                         SET FLTBEG TO FLTBEG-(FLTCNT OF LAST BLOCK)
*                     ENDIF
*                     IF FIRST FILTERING FOR BLOCK AND (PRESENT BLOCK
*                       IS SMALLER THAN (NUMCOF-1)/2)
*                     THEN
*                         SET FLTBEG TO FLTBEG-( (NUMCOF-1)/2)
*                     ENDIF
*                     CALL FILTDAT TO FILTER TRIPLETS SINCE LAST CALL UP TO
*                       PRESENT FILL GAP
*                     IF JCONT EQUALS 0
*                     THEN
*                         INCREMENT ISTATS(4) BY 1
*                     ENDIF
*                     SET INTCNT TO 0
*                     SET FLTBEG TO COUNTR
*                     SET FLTCNT TO 1
*                 ELSE
*                     SET LARGE FLAG TO 0
*                     CALL INTERP TO INTERPOLATE
*                     IF LARGE FLAG EQUAL 1
*                     THEN
*                         INCREMENT INTCNT BY 1
*                     ELSE
*                         SET INTCNT TO 0
*                     ENDIF
*                 ENDIF
*             ENDIF
*             INCREMENT COUNTR BY 1
*         ENDDO
*         RETURN
*******************************************************************
      REAL*4 TRIPS(3,0:19999), FTRIPS(3,0:19999), COEF(500)
      INTEGER*4 ISTART,INUM,ISTATS(10),NUMCOF,DATGAP,INTPTS,ICONT
      INTEGER*4 FLTBEG,FLTCNT,COUNTR,INTCNT,INTBEG,JCONT,BCONT,FITDEG
      INTEGER*4 OVRFLT,OVERRN/0/,UNDRRN/0/
C
      DATA FILL/999.0/,IFIRST/0/
C
      IF(IFIRST.EQ.0) THEN
         IFIRST = 1
         BCONT = 0
         UNDRRN = 0
      ENDIF
C
      OVRFLT = 0
      FLTBEG = ISTART
      FLTCNT = 0
C
      COUNTR = ISTART
      INTCNT = 0
C
      IAGAIN = 0
C
      PRINT *, ' IN FINDFILL:'
      PRINT *,' NUMCOF,FLTBEG,INUM,DATGAP,INTPTS,ISTATS,ICONT'
      PRINT *,NUMCOF,FLTBEG,INUM,DATGAP,INTPTS
      PRINT *,ISTATS,ICONT
C
C     REPEAT FOR THE NUMBER OF TRIPLETS IN THIS GROUP
C
      DO 5 I=1,INUM
C
C     INITIATE OUTPUT TRIPLETS TO FILL
C
          DO 6 K=1,3
             FTRIPS(K,MOD(COUNTR,20000)) = FILL
   6      CONTINUE
C
          FLTCNT = FLTCNT + 1
C
C     CHECK FOR FILL DATA OR END OF DATA...KEEP COUNT
C
          IF (TRIPS(1,MOD(COUNTR,20000)).EQ.FILL.OR.
     I        TRIPS(2,MOD(COUNTR,20000)).EQ.FILL.OR.
     I        TRIPS(3,MOD(COUNTR,20000)).EQ.FILL.OR.I.EQ.INUM)THEN
C
C     SET CONTINUATION FLAG
C
             IF(I.EQ.INUM.AND.ICONT.EQ.1) THEN
                 JCONT = 1
             ELSE
                 JCONT = 0
             ENDIF
C
             IF (INTCNT.EQ.0) THEN
                 INTBEG = COUNTR
             ENDIF
C
             INTCNT = INTCNT + 1
          ENDIF
C
C     CHECK SITUATION OF FILL DATA AT END OF BLOCK OF TRIPLETS WHICH
C     FOLLOWED A DATA GAP...GO FILTER DATA.
C
          IF(JCONT.EQ.0.AND.I.EQ.INUM) THEN
              IF(TRIPS(1,MOD(COUNTR,20000)).EQ.FILL.OR.
     I          TRIPS(2,MOD(COUNTR,20000)).EQ.FILL.OR.
     I          TRIPS(3,MOD(COUNTR,20000)).EQ.FILL) THEN
C
                    IF (IAGAIN.EQ.1) FLTCNT = FLTCNT + 1
C
                    FLTCNT = FLTCNT - INTCNT
      PRINT *,' IN JCONT=0,I=INUM: FLTBEG,FLTCNT,JCONT,ISTATS,'
      PRINT *,' INTCNT,I:'
      PRINT *,FLTBEG,FLTCNT,JCONT,ISTATS,INTCNT,I
C
CCCCC ADDED FOR VERSION 2.1 - TO CHECK BEGIN OF DATA WHEN FILTERING IS
CCCCC DUE TO FILL AT END OF BLOCK.
C
                    OVERRN = FLTCNT - (NUMCOF-1)/2
C
C     ADJUST NUMBER OF TRIPLETS TO FILTER IF THERE ARE LESS TRIPLETS
C     IN THE BLOCK BEFORE THE NEXT BREAK THAN THERE ARE 1/2 THE
C     NUMBER OF COEFFICIENTS.
C
                    IF(FLTBEG.EQ.ISTART.AND.UNDRRN.LT.0.AND.
     I                      BCONT.EQ.1) THEN
                         FLTBEG = FLTBEG + UNDRRN
                         FLTCNT = FLTCNT - UNDRRN
                         BCONT=0
                         IF(FLTBEG.EQ.ISTART.AND.OVERRN.LT.0.AND.
     I                         BCONT.EQ.1) THEN
                             OVERRN = 0
                         ENDIF
                     ELSE
C
CCCCCCC ADDED IN V.2.2 - BACKUP IF SMALL AMOUNT OF FILL FOUND AFTER A BREAK
C
                         IF (OVRFLT.LT.0) THEN
                             FLTBEG = FLTBEG + OVRFLT
                             FLTCNT = FLTCNT - OVRFLT
                             IF(FLTBEG.EQ.ISTART.AND.OVERRN.LT.0.AND.
     I                             BCONT.EQ.1) THEN
                                 OVERRN = 0
                             ENDIF
                         ELSE
CCCCCCC
C
C     ADJUST AMOUNT OF DATA TO BE FILTERED IF THERE ARE FEWER TRIPLETS
C     THAN THERE ARE 1/2 NUMBER OF COEFFICIENTS, AT BEGIN OF A BLOCK
C     AND BEFORE TIME GAP.
C
                             IF(FLTBEG.EQ.ISTART.AND.OVERRN.LT.0.AND.
     I                             BCONT.EQ.1) THEN
                                FLTBEG = FLTBEG + OVERRN
                                FLTCNT = FLTCNT - OVERRN
                                OVERRN = 0
                             ENDIF
                         ENDIF
                     ENDIF
                     OVRFLT = 0
C
C
      PRINT *,' BEFORE CALLING FILTDAT: FLTBEG,FLTCNT,JCONT,ISTATS,'
      PRINT *,' INTCNT,I,ISTART:'
      PRINT *,FLTBEG,FLTCNT,JCONT,ISTATS,INTCNT,I,ISTART
C
                    CALL FILTDAT (TRIPS, FTRIPS, COEF, NUMCOF,
     C                   FLTBEG, FLTCNT, JCONT, BCONT, ISTATS)
C
C
                    IF(I.EQ.INUM.AND.OVERRN.LT.0.AND.JCONT.EQ.1) THEN
                        UNDRRN = -FLTCNT
                    ELSE
                        UNDRRN = 0
                    ENDIF
C
CCCCC END OF CODE ADDED FOR VERSION 2.1
CCCCC
C
                    IF (JCONT.EQ.0) THEN
                        ISTATS(4) = ISTATS(4) + 1
                    ENDIF
C
                    INTCNT = 0
                    FLTCNT = 1
                    IAGAIN = 1
                    FLTBEG = COUNTR
              ENDIF
          ENDIF
C
C     INTCNT GT 0 MEANS FILL ENCOUNTERED...DO THIS WHEN FIRST GOOD TRIPLET
C     IS ENCOUNTERED.
C
          IF (INTCNT.GT.0.AND.TRIPS(1,MOD(COUNTR,20000)).NE.FILL.AND.
     I        TRIPS(2,MOD(COUNTR,20000)).NE.FILL.
     I        AND.TRIPS(3,MOD(COUNTR,20000)).NE.FILL) THEN
C
      PRINT *,' I,INUM,FLTCNT,COUNTR,INTBEG,INTCNT,TRIPS(COUNTR)'
      PRINT *,I,INUM,FLTCNT,COUNTR,INTBEG,INTCNT,
     P        TRIPS(1,MOD(COUNTR,20000))
C
C     CALL FILTDAT DUE TO LARGE FILL DATA GAP OR END OF BLOCK
C
             IF (INTCNT+2.GT.DATGAP. OR .I.EQ.INUM) THEN
C
                 IF (I.NE.INUM) THEN
                     FLTCNT = FLTCNT - 1
                 ENDIF
C
                 IF (I.EQ.INUM)
     I                  FLTCNT = FLTCNT + 1
C
                 FLTCNT = FLTCNT - INTCNT
C
                 OVERRN = FLTCNT - (NUMCOF-1)/2
C
C     ADJUST NUMBER OF TRIPLETS TO FILTER IF THERE ARE LESS TRIPLETS
C     IN THE BLOCK BEFORE THE NEXT BREAK THAN THERE ARE 1/2 THE
C     NUMBER OF COEFFICIENTS.
C
                 IF(FLTBEG.EQ.ISTART.AND.UNDRRN.LT.0.AND.
     I                  BCONT.EQ.1) THEN
                     FLTBEG = FLTBEG + UNDRRN
                     FLTCNT = FLTCNT - UNDRRN
                     BCONT=0
                     IF(FLTBEG.EQ.ISTART.AND.OVERRN.LT.0.AND.
     I                     BCONT.EQ.1) THEN
                         OVERRN = 0
                     ENDIF
                 ELSE
C
CCCCCCC ADDED IN V.2.2 - BACKUP IF SMALL AMOUNT OF FILL FOUND AFTER A BREAK
C
                     IF (OVRFLT.LT.0) THEN
                         FLTBEG = FLTBEG + OVRFLT
                         FLTCNT = FLTCNT - OVRFLT
                         IF(FLTBEG.EQ.ISTART.AND.OVERRN.LT.0.AND.
     I                         BCONT.EQ.1) THEN
                             OVERRN = 0
                         ENDIF
                     ELSE
CCCCCCC
C
C     ADJUST AMOUNT OF DATA TO BE FILTERED IF THERE ARE FEWER TRIPLETS
C     THAN THERE ARE 1/2 NUMBER OF COEFFICIENTS, AT BEGIN OF A BLOCK
C     AND BEFORE TIME GAP.
C
                         IF(FLTBEG.EQ.ISTART.AND.OVERRN.LT.0.AND.
     I                         BCONT.EQ.1) THEN
                            FLTBEG = FLTBEG + OVERRN
                            FLTCNT = FLTCNT - OVERRN
                            OVERRN = 0
                         ENDIF
                     ENDIF
                 ENDIF
                 OVRFLT = 0
C
C
      PRINT *,' BEFORE CALLING FILTDAT: FLTBEG,FLTCNT,JCONT,ISTATS,'
      PRINT *,' INTCNT,I,ISTART:'
      PRINT *,FLTBEG,FLTCNT,JCONT,ISTATS,INTCNT,I,ISTART
C
                 CALL FILTDAT (TRIPS, FTRIPS, COEF, NUMCOF,
     C                FLTBEG, FLTCNT, JCONT, BCONT, ISTATS)
C
C
                 IF(I.EQ.INUM.AND.OVERRN.LT.0.AND.JCONT.EQ.1) THEN
                     UNDRRN = -FLTCNT
                 ELSE
                     UNDRRN = 0
                 ENDIF
C
C
                 IF (JCONT.EQ.0) THEN
                     ISTATS(4) = ISTATS(4) + 1
                 ENDIF
C
                 INTCNT = 0
                 FLTCNT = 1
                 IAGAIN = 1
                 FLTBEG = COUNTR
C
             ELSE
C
C     ADJUST INTERP COUNTERS TO INCLUDE LAST AND FIRST GOOD POINTS (AS
C     LONG AS THEY ARE NOT BEYOND ENDS OF BLOCK)
C
                 INTBEG = INTBEG - 1
                 INTCNT = INTCNT + 2
C
C     CALL INTERP TO INTERPOLATE SMALL DATA GAPS
C
                 LARGE = 0
C
      PRINT *,' BEFORE INTERP: I,INUM,INTCNT,INTBEG,INTPTS,LARGE,
     P                         FLTBEG,FLTCNT:'
      PRINT *, I,INUM,INTCNT,INTBEG,INTPTS,LARGE,FLTBEG,FLTCNT
C
                 IF (INTBEG.LE.FLTBEG+INTPTS.AND.BCONT.EQ.0) THEN
C
                     FLTBEG = INTBEG + INTCNT
                     FLTCNT = 0
                 ELSE
C
CCCCCC ADDED IN V.2.2 - TO DETERMINE SMALL AMT OF FILL JUST AFTER BREAK
C
                     IF (FLTBEG.EQ.ISTART.AND.BCONT.EQ.1.AND.
     I                           OVRFLT.EQ.0) THEN
                        OVRFLT = FLTCNT - (NUMCOF-1)/2 - INTCNT - 1
                     ENDIF
CCCCCCC
                     CALL INTERP (TRIPS, INTCNT, INTBEG, INTPTS, LARGE,
     C                DATGAP, ISTATS, FITDEG)
C
                 ENDIF
C
      PRINT *,' AFTER INTERP: I,INUM,INTCNT,INTBEG,INTPTS,LARGE,
     P          FLTBEG,FLTCNT,ISTART,OVRFLT:'
      PRINT *, I,INUM,INTCNT,INTBEG,INTPTS,LARGE,FLTBEG,FLTCNT,
     P          ISTART,OVRFLT
C
                 IF (LARGE.EQ.1) THEN
                     INTCNT = INTCNT + 1
                 ELSE
                     INTCNT = 0
                 ENDIF
             ENDIF
          ENDIF
          COUNTR = COUNTR +1
  5   CONTINUE
      RETURN
      END
C
C
C
C
C
C
      SUBROUTINE INTERP (TRIPS, JNTCNT,INTBEG, INTPTS,
     S                   LARGE, DATGAP, ISTATS, FITDEG)
*************************************************************************
*
*      TITLE:  INTERP FILL DATA
*
*      FILE NAME:  INTERP.FOR
*
*      PURPOSE:  EXAMINES TRIPLETS ADJACENT TO "SMALL FILL DATA GAP" AND
*                         SETS DATA UP FOR SPLINE INTERPOLATION.
*
*      HISTORY:
*
*        AUTHOR               DATE            CHANGE
*        ------               ----            ------
*        S.J.KEMPLER          12/1/88         ORIGINAL CODE
*                             2/16/89         IMPLEMENT LEAST SQUARES INTER-
*                                             POLATION AND FITDEG INPUT.
*
*      CALLING SEQUENCE:  SUBROUTINE INTERP (TRIPS, JNTCNT,INTBEG, INTPTS,
*                             LARGE, DATGAP, ISTATS, FITDEG)
*
*      ARGUMENT LIST:
*
*       NAME               TYPE   USE    PURPOSE
*       ----               ----   ---    -------
*
*	TRIPS(3,0:19999)   R*4     I     DATA TO BE FILTERED
*       JNTCNT             I*4     I     NUMBER OF TRIPLETS TO BE INTERPOLATED
*                                        OVER
*       INTBEG             I*4     I     NUMBER OF FILL TRIPLET TO START
*                                        INTERPOLATION
*       INTPTS             I*4     I     THE NUMBER OF ADJACENT TRIPLETS USED IN
*                                        THE INTERPOLATION
*       LARGE              I*4     O     FLAG PASSED BACK INDICATING THAT THE
*                                        "SMALL" AMOUNT OF FILL DATA HAS BECOME
*                                        "LARGE" DUE TO INSUFFICIENT NUMBER OF
*                                        GOOD DATA ADJACENT TO THE FILL DATA,
*                                        THEREBY FORCING THE GOOD ADJACENT DATA
*                                        TO BE MADE FILL.  0: SMALL; 1: MADE
*                                        LARGE.
*       DATGAP             I*4     I     THE NUMBER OF CONSECUTIVE FILL DATA
*                                        VALUES (999.0) PERMITTED BEFORE THE
*                                        "GAP" IS TOO LARGE TO INTERPOLATE.
*                                        (24 MAX)
*       ISTATS(10)         I*4     O     ISTATS(1)=NUMBER OF INTERPOLATED
*                                        TRIPLETS ADDED; ISTATS(2)=NUMBER OF
*                                        INTERPOLATED PERIODS; ISTATS(3)=N0.OF
*                                        FILL TRIPLETS ADDED; ISTATS(4)= NUMBER
*                                        OF GAPS AND LARGE MISSING DATA
*                                        PERIODS.
*                                        ISTATS(5)=TOTAL NUMBER OF TRIPS.
*        FITDEG            I*4     I     DEGREE OF LEAST SQUARES FIT
*
*      MODULES REFERENCED:  LSPOL
*
*      NOTES:  EXAMINES TRIPLETS 'INTPTS' TRIPLETS AFTER THE FILL TO
*             ASSURE THAT THEY ARE GOOD.  IF NOT, INCREASE FILL GAP TO INCLUDE
*             FILL DATA AND THE GOOD DATA IN BETWEEN.  THIS MAY MAKE THE FILL
*             GAP LARGE> SET LARGE FLAG.  IF NOT, SET UP AND SPLINE INTERPOLATE.
*
*      PDL:
*
*         SET INTCNT EQUAL TO JNTCNT
*         SET ITER TO -1
*         DO WHILE ITER IS LESS THAN INTPTS
*             INCREMENT ITER BY 1
*             IF INTCNT IS GREATER THAN DATGAP
*             THEN
*                 SET LARGE FLAG TO 1
*                 RETURN
*             ELSE
*                 IF TRIPS(INTBEG+INTCNT+ITER) IS EQUAL TO FILL
*                 THEN
*                     DOFOR M EQUAL TO (INTBEG+INTCNT) TO (INTBEG+INTCNT+
*                      ITER-1)
*                         SET TRIPS(M) TO FILL
*                     ENDDO
*                     SET INTCNT EQUAL TO INTCNT + ITER
*                     SET ITER TO -1
*                 ENDIF
*             ENDIF
*         ENDDO
*         SET ICT TO 0
*         DO WHILE ICT IS LESS THAN INTPTS
*             SET TRIPS(INTBEG-INTPTS+ICT) TO ADJTRP(ICT+1)
*             SET (INTBEG-INTPTS+ICT) TO ADJTIM(ICT+1)
*             SET TRIPS(INTBEG+INTCNT+ICT) TO ADJTRP(INTPTS+ICT+1)
*             SET (INTBEG+INTCNT+ICT) TO ADJTIM(INTPTS+ICT+1)
*             INCREMENT ICT BY 1
*         ENDDO
*         CALL LSPOL TO INITIALIZE THE SPLINE INTERPOLATION
*         INCREMENT ISTATS(2) BY 1
*         SET ICT TO 0
*         DO WHILE ICT IS LESS THAN INTCNT
*             SET INTBEG + ICT TO NEXTPT
*             INTERPOLATE FOR CURRENT TRIPLET NUMBER
*             INCREMENT ISTATS(1) BY 1
*         ENDDO
*         RETURN
*******************************************************************
      REAL*4 TRIPS(3,0:19999)
      REAL ADJTIM(100), ADJTRP(100), ERR(100)
      INTEGER*4 ISTATS(10),INTPTS,DATGAP,INTCNT,INTBEG,ITER,JNTCNT
      INTEGER*4 FITDEG
      DIMENSION CCC(6), AAA(1800), FFF(100), KKK(24)
C
      DATA FILL/999.0/
C
      INTCNT = JNTCNT
C
      NON = 0
      ICOUNT = 0
      ITER=0
C
 12   IF (ITER.GE.INTPTS) GO TO 15
C
          ICOUNT = ICOUNT + 1
C
C     CHECK FOR ENLARGED DATA GAP
C
          IF (INTCNT.GT.DATGAP) THEN
              LARGE = 1
              RETURN
          ELSE
C
C     EXAMINE TRIPLETS ON EITHER SIDE OF SMALL GAP FOR FILL.  IF FOUND
C     MAKE ALL THE TRIPLETS UP TO THAT FILLED TRIPLET FILL AS WELL.  GO
C     BACK AND RE-EXAMINE THE SIZE UP THE ENLARGED DATA GAP.  IF STILL
C     SMALL, REPEATTHIS ACTION.
C
              INTRP = INTBEG + INTCNT + ITER + NON
              IF (TRIPS(1,MOD(INTRP,20000)).EQ.FILL
     I                 .OR.TRIPS(2,MOD(INTRP,20000)).EQ.FILL
     I                 .OR.TRIPS(3,MOD(INTRP,20000)).EQ.FILL)   THEN
                  IF(ITER.GT.0) THEN
                 DO 5 M = INTBEG+INTCNT, INTRP-1
                         DO 10 J=1,3
                             TRIPS(J,MOD(M,20000)) = FILL
   10                    CONTINUE
    5                CONTINUE
                  ENDIF
C
                  INTCNT = INTCNT + ITER + NON
                  ITER = -1
                  NON = 1
              ENDIF
          ENDIF
          ITER = ITER + 1
C      PRINT *,' INTERP: INTBEG,INTCNT,ITER,INTRP,INTPTS,ICOUNT,FITDEG:'
C      PRINT *, INTBEG,INTCNT,ITER,INTRP,INTPTS,ICOUNT,FITDEG
      GO TO 12
C
C     FINALLY, AFTER ALL THE REQUIRED TRIPLETS ARE EXAMINED ON EITHER
C     OF THE GAP, SET UP FOR THE CALL TO THE SPLINE INTERPOLATION, FOR
C     EACH AXIS.
C
  15  DO 20 J=1,3
C
      ICT = 1
C
      DO 16 K=1,INTPTS
C
          ADJTRP(ICT) = TRIPS(J,MOD(INTBEG-INTPTS+ICT-1,20000))
          ADJTIM(ICT) = -(INTPTS)+ICT-1
C
      ADJTRP(INTPTS+ICT) = TRIPS(J,MOD(INTBEG+INTCNT+ICT+NON-1,20000))
      ADJTIM(INTPTS+ICT) = INTCNT+ICT+NON-1
C
          ICT = ICT + 1
C
C      PRINT *,' ICT,(ADJTIM(KK),ADJTRP(KK),KK=1,2*INTPTS),
C     P          INTPTS,DS1,DSM:'
C      PRINT *,ICT,(ADJTIM(KK),ADJTRP(KK),KK=1,2*INTPTS),
C     P          INTPTS,DS1,DSM
C
  16  CONTINUE
C
C      DS1 = 1.0E30
C      DSM = 1.0E30
C
      PRINT *,' J,ICT,(ADJTIM(KK),ADJTRP(KK),KK=1,2*INTPTS),
     P    INTPTS,DS1,DSM:'
      PRINT *,J,ICT,(ADJTIM(KK),ADJTRP(KK),KK=1,2*INTPTS),
     P          INTPTS,DS1,DSM
C
C     CALL R. THOMPSON'S SPLINE_INIT AND SPLINE.
C      CALL SPLINE_INIT (ADJTIM,ADJTRP,2*INTPTS,DS1,DSM,DDS)
C
C     VERSION 2.0: CALL R. THOMPSON'S LSPOL ROUTINE INSTEAD
C
      CALL LSPOL (ADJTRP,ADJTIM,2*INTPTS,FITDEG,CCC,RMS,ERR,AAA,FFF,KKK)
C
      IF(J.EQ.1) ISTATS(2) = ISTATS(2) + 1
      ICT = 0
C
      DO 18 K=1,INTCNT+NON
          NEXTPT = INTBEG + ICT
C
C     VERSION 2.0: SPLINE FUNCTION COMMENTED OUT, ITERATION ADDED
C
C          TRIPS(J,MOD(INTBEG+ICT,20000)) =
C     C          SPLINE(ADJTIM,ADJTRP,2*INTPTS,DDS,FLOATJ(INTBEG+ICT))
C
      TRIPS(J,MOD(NEXTPT,20000)) = CCC(1)
C
          DO 17 KC=2,FITDEG+1
C
             TRIPS(J,MOD(NEXTPT,20000)) = TRIPS(J,MOD(NEXTPT,20000)) +
     T          ( CCC(KC) * ( FLOAT(ICT)**(KC-1) ) )
C
  17      CONTINUE
          IF(J.EQ.1) THEN
             ISTATS(1) = ISTATS(1) + 1
             PRINT *,' TRIPS(J,MOD(NEXTPT,20000)):',
     P                 TRIPS(J,MOD(NEXTPT,20000))
          ENDIF
          ICT = ICT +1
  18  CONTINUE
C
   20 CONTINUE
C
C      PRINT *,' ICOUNT:',ICOUNT
C
      RETURN
      END
C
C
C
C
      SUBROUTINE FILTDAT (TRIPS,FTRIPS,COEF,NUMCOF,FLTBEG,
     S                FLTCNT,JCONT,BCONT,ISTATS)
*************************************************************************
*
*      TITLE:  FILTER DATA
*
*      FILE NAME:  FILTDAT.FOR
*
*      PURPOSE:  RECIEVES THE DATA TO BE FILTERED AND RETURNS THE FILTERED
*                  DATA
*
*      HISTORY:
*
*        AUTHOR               DATE            CHANGE
*        ------               ----            ------
*        S.J.KEMPLER          12/1/88         ORIGINAL CODE
*
*      CALLING SEQUENCE:  SUBROUTINE FILTDAT (TRIPS, FTRIPS, COEF, NUMCOF,
*                           FLTBEG, FLTCNT, JCONT, BCONT, ISTATS)
*
*      ARGUMENT LIST:
*
*       NAME             TYPE   USE    PURPOSE
*       ----             ----   ---    -------
*
*	TRIPS(3,0:19999)   R*4     I     DATA TO BE FILTERED.
*	FTRIPS(3,0:19999)  R*4     O     FILTERED DATA.
*       COEF(500)          R*4     O     THE ENTIRE ARRAY OF FILTER COEF TO BE
*                                        USED. (2N-1, IF SYMMETRICAL; N IF NOT
*                                        SYMMETRICAL, WHERE N IS THE NUMBER
*                                        OF COEFF RECIEVED FROM THE FILTER COEF
*                                        FILE.  ONLY HALF OF THE COEF ARE PRO-
*                                        VIDED IF THE FILTER IS SYMMETRICAL.)
*       NUMCOF             I*4     I     TOTAL NUMBER OF FILTER COEF (SEE
*                                        ABOVE).
*       FLTBEG             I*4     I     NUMBER OF FIRST TRIPLET TO BE
*                                        FILTERED.
*       FLTCNT             I*4     I     NUMBER OF TRIPLETS TO BE FILTERED.
*       JCONT              I*4     I     CONTINUATION FLAG FOR END OF
*                                        DATA GROUP: 0=NEXT GROUP
*                                        OF TRIPLETS FOLLOWS DATA GAP; 1=
*                                        NEXT GROUP OF TRIPLETS IS A
*                                        CONTINUATION OF GOOD DATA FROM
*                                        LAST GROUP.
*       BCONT              I*4     I     CONTINUATION FLAG FOR BEGINNING
*                                        OF NEXT DATA GROUP: 0=NEXT GROUP
*                                        OF TRIPLETS FOLLOWS DATA GAP; 1=
*                                        NEXT GROUP OF TRIPLETS IS A
*                                        CONTINUATION OF GOOD DATA FROM
*                                        LAST GROUP.
*        ISTATS(10)        I*4   I/O     ISTATS(1)=NUMBER OF INTERPOLATED
*                                        TRIPLETS ADDED; ISTATS(2)=NUMBER OF
*                                        INTERPOLATED PERIODS; ISTATS(3)=N0.OF
*                                        FILL TRIPLETS ADDED; ISTATS(4)= NUMBER
*                                        OF GAPS AND LARGE MISSING DATA PERIODS.
*                                        ISTATS(5)=TOTAL NUMBER OF TRIPS.
*
*      MODULES REFERENCED:  FILDAT, FILTPT
*
*      NOTES:  FIRST, BE SURE THERE ARE MORE POINTS TO BE FILTERED THAN
*              THERE ARE COEFFICIENTS TO DO THE FILTERING.  THEN, IF
*              THIS IS NOT A CONTINUATION OF REPLACE ONE
*              HALF THE NUMBER OF COEFFICIENTS OF DATA AT THE BEGINNING
*              AND END OF THIS BLOCK OF DATA.  INBETWEEN, CALL FILTPT
*
*
*      PDL:
*
*         IF IFIRST EQUALS 0
*         THEN
*            SET BCONT EQUAL TO 0 (DATA NOT BEING CONTINUED)
*            SET IFIRST EQUAL TO 1
*         ENDIF
*         IF FLTCNT IS LESS THAN NUMCOF AND BCONT EQUAL 1 AND
*           JCONT EQUALS 1
*         THEN
*            CALL FILDAT TO FILL FLTCNT NUMBER OF TRIPLETS
*            RETURN
*         ENDIF
*         SET IBEGIN EQUAL TO FLTBEG
*         SET ICOUNT EQUAL TO FLTCNT
*         IF BCONT EQUALS 0
*         THEN
*            CALL FILDAT TO FILL THE FIRST (NUMCOF-1)/2 WITH FILL
*            SET IBEGIN EQUAL TO IBEGIN + (NUMCOF-1)/2
*            SET ICOUNT EQUAL TO ICOUNT - (NUMCNT-1)/2
*         ENDIF
*         IF JCONT EQUALS 0
*         THEN
*            CALCULATE THE BEGINNING OF THE LAST (NUMCOF-1)/2
*            CALL FILDAT TO FILL THE LAST (NUMCOF-1)/2 WITH FILL
*            SET ICOUNT EQUAL TO ICOUNT - (NUMCNT-1)/2
*         ENDIF
*         CALL FILTPT TO FILTER ACROSS DATA THAT HAS A SUFFICIENT
*           NUMBER OF ADJACENT TRIPLETS
*         SET JCONT TO BCONT
*         RETURN
*******************************************************************
      REAL*4 TRIPS(3,0:19999),FTRIPS(3,0:19999),COEF(500)
      INTEGER*4 NUMCOF,FLTBEG,FLTCNT,JCONT,BCONT,ISTATS(10)
C
      DATA IFIRST/0/
C
      PRINT *,' IN FILTDAT:'
C
      IF (IFIRST.EQ.0) THEN
          IFIRST = 1
          BCONT = 0
      ENDIF
C
C     IF THE AMOUNT OF DATA IS LESS THAN THE NUMBER OF COEFFICIENTS,
C     MAKE IT ALL FILL.
C
      IF(FLTCNT.LT.NUMCOF.AND.BCONT.EQ.0.AND.JCONT.EQ.0) THEN
C
      PRINT *,' FLTCNT<NUMCOF: FLTCNT,NUMCOF:',FLTCNT,NUMCOF
C
          CALL FILDAT (FTRIPS,FLTBEG,FLTCNT,ISTATS)
          RETURN
C
      ENDIF
C
      IBEGIN = FLTBEG
      ICOUNT = FLTCNT
C
C     FILL THE BEGINNING OF THE SET OF TRIPLETS IF FOLLOWING A GAP
C
      IF (BCONT.EQ.0) THEN
C
      PRINT *,' ICOUNT,IBEGIN,(NUMCOF-1)/2,NUMCOF:'
      PRINT *,ICOUNT,IBEGIN,(NUMCOF-1)/2,NUMCOF
C
          CALL FILDAT(FTRIPS,IBEGIN,(NUMCOF-1)/2,ISTATS)
          IBEGIN = IBEGIN + (NUMCOF-1)/2
          ICOUNT = ICOUNT - (NUMCOF-1)/2
      ENDIF
C
C     FILL THE END OF THE SET OF TRIPLETS IF A GAP FOLLOWS
C
      IF (JCONT.EQ.0) THEN
          IBGEND = FLTBEG+FLTCNT-(NUMCOF-1)/2
          CALL FILDAT(FTRIPS,IBGEND,(NUMCOF-1)/2,ISTATS)
          ICOUNT = ICOUNT - (NUMCOF-1)/2
      ENDIF
C
      PRINT *,' BEFORE FILTPT: ICOUNT,IBEGIN,(NUMCOF-1)/2,NUMCOF,
     P          BCONT,JCONT:'
      PRINT *,ICOUNT,IBEGIN,(NUMCOF-1)/2,NUMCOF,BCONT,JCONT
C
C     FINALLY...
C
      CALL FILTPT (TRIPS,FTRIPS,IBEGIN,ICOUNT,COEF,NUMCOF)
      BCONT = JCONT
C
      RETURN
      END
C
C
C
C
C
      SUBROUTINE FILDAT(FTRIPS,ITRIP1,NTRIPS,ISTATS)
C
*************************************************************************
*
*      TITLE:  ENTER FILL DATA
*
*      FILE NAME:  FILDAT.FOR
*
*      PURPOSE:  LOADS FILL DATA INTO OUT FILTERED DATA BUFFER
*                ACCORDING TO FIRST TRIPLET AND NUMBER OF TRIPLETS.
*
*      HISTORY:
*
*        AUTHOR               DATE            CHANGE
*        ------               ----            ------
*        S.J.KEMPLER          12/1/88         ORIGINAL CODE
*
*      CALLING SEQUENCE: SUBROUTINE FILDAT (FTRIPS, ITRIP1, NTRIPS, ISTATS)
*
*      ARGUMENT LIST:
*
*       NAME             TYPE   USE    PURPOSE
*       ----             ----   ---    -------
*
*       FTRIPS(3,0:19999) R*4   I/O    FILTERED DATA BUFFER WHERE FILL
*                                      DATA IS PLACED AS DESIRED.
*       ITRIP1           I*4     I     THE NUMBER OF THE FIRST TRIPLET
*                                      TO RECEIVE FILL.
*       NTRIPS           I*4     I     THE NUMBER OF TRIPLETS TO
*                                      RECEIVE FILL.
*      ISTATS(10)        I*4    I/O     ISTATS(1)=NUMBER OF INTERPOLATED
*                                      TRIPLETS ADDED; ISTATS(2)=NUMBER OF
*                                      INTERPOLATED PERIODS; ISTATS(3)=N0.OF
*                                      FILL TRIPLETS ADDED; ISTATS(4)= NUMBER
*                                      OF GAPS AND LARGE MISSING DATA PERIODS.
*                                      ISTATS(5)=TOTAL NUMBER OF TRIPS.
*
*      MODULES REFERENCED:  NONE
*
*
*      NOTES:
*
*      THIS ROUTINE ISOLATES ALL OF THE TRIPLET FILLING ACTION
*
*      PDL:
*
*         SET COUNTR TO 0
*         DO WHILE COUNTR IS LESS THAN OR EQUAL TO NTRIPS
*            SET FTRIPS(ITRIP1+COUNTR) TO FILL
*            INCREMENT COUNTR BY 1
*         ENDDO
*         RETURN
*******************************************************************
C
      REAL FTRIPS(3,0:19999)
      INTEGER*4 ITRIP1,NTRIPS,ISTATS(10)
C
      DATA FILL/999.0/
C
      DO 10 J=ITRIP1,ITRIP1+NTRIPS-1
          DO 5 I=1,3
               FTRIPS(I,MOD(J,20000)) = FILL
   5      CONTINUE
          ISTATS(3) = ISTATS(3) + 1
  10  CONTINUE
C
      RETURN
      END
C
C
C
C
      SUBROUTINE FILTPT (TRIPS,FTRIPS,ITRIP1,NTRIPS,COEF,NUMCOF)
C
*************************************************************************
*
*      TITLE:  FILTER POINTS
*
*      FILE NAME:  FILTPT.FOR
*
*      PURPOSE:  PERFORMS FILTERING ON THE GIVEN CONSECUTIVE TRIP[LETS
*
*      HISTORY:
*
*        AUTHOR               DATE            CHANGE
*        ------               ----            ------
*        S.J.KEMPLER          12/1/88         ORIGINAL CODE
*
*      CALLING SEQUENCE: SUBROUTINE FILTPT (TRIPS, FTRIPS, IBEGIN, ICOUNT,
*                           COEF, NUMCOF)
*
*      ARGUMENT LIST:
*
*       NAME             TYPE   USE    PURPOSE
*       ----             ----   ---    -------
*
*       TRIPS(3,0:19999) R*4     I     DATA TO BE FILTERED
*       FTRIPS(3,0:19999) R*4   I/O    BUFFER TO RECEIVE FILTERED DATA
*       ITRIP1           I*4     I     THE NUMBER OF THE FIRST TRIPLET
*                                      TO BE FILTERED
*       NTRIPS           I*4     I     THE NUMBER OF TRIPLETS TO BE
*                                      FILTERED
*       COEF(500)        R*4     O     THE ENTIRE ARRAY OF FILTER COEF TO BE
*                                      USED. (2N-1, IF SYMMETRICAL; N IF NOT
*                                      SYMMETRICAL, WHERE N IS THE NUMBER
*                                      OF COEFF RECIEVED FROM THE FILTER COEF
*                                      FILE.  ONLY HALF OF THE COEF ARE PRO-
*                                      VIDED IF THE FILTER IS SYMMETRICAL.)
*       NUMCOF           I*4     I     TOTAL NUMBER OF FILTER COEF (SEE
*                                      ABOVE).
*
*
*      MODULES REFERENCED:  NONE
*
*
*      NOTES:
*
*      ALL OF THE TRIPLET LOCATIONS TO BE FILTERED FOR ARE DONE SO HERE.
*      THE BEGINNING AND END OF DATA SEGMENTS, ALTHOUGH NATURALLY USED
*      FOR THE FILTERING, HAVE FILL IN THE CORRESPONDING OUTPUT
*      LOCATIONS, SINCE THOSE TRIPLETS CAN'T BRE FILTERED.  THIS FILLING
*      HAS ALREADY BEEN DONE IN SUBROUTINE FILDAT.
*
*      PDL:
*
*         DOFOR I EQUAL TO 1 TO NTRIPS
*             INITIALIZE FTRIPS(ITRIP1+I-1) TO 0
*             DOFOR J EQUAL TO 1 TO NUMCOF
*                 SET FTRIPS(ITRIP1+I-1) TO FTRIPS(ITRIP1+I-1) +
*                   COEF(J) * TRIPS(ITRIP1 + J - (NUMCOF-1)/2 -1)
*             ENDDO
*         ENDDO
*         RETURN
*******************************************************************
C
      REAL*4 TRIPS(3,0:19999),FTRIPS(3,0:19999),COEF(500)
      INTEGER*4 ITRIP1,NTRIPS,NUMCOF
C
      NCFHAF = (NUMCOF-1)/2
C
      PRINT *,' IN FILTPT: NCFHAF:',NCFHAF
C
C      PRINT *,' FTRIPS BEFORE:',(JJ,FTRIPS(1,JJ),JJ=0,19999,20)
C
C      FOR EACH TRIPLET...
C
      DO 5 I=1,NTRIPS
          II = I - 1
C
C      FOR EACH AXIS...
C
          DO 10 J=1,3
              FTRIPS(J,MOD(ITRIP1+II,20000)) = 0.0
C
C      FOR EACH COEFFICIENT...
C
              DO 15 K=1,NUMCOF
                  FTRIPS(J,MOD(ITRIP1+II,20000)) =
     F            FTRIPS(J,MOD(ITRIP1+II,20000)) +
     F        (TRIPS(J,MOD(ITRIP1+II + K-1 - NCFHAF,20000)) * COEF(K))
  15          CONTINUE
C
C      PRINT *,' J,I,TRIPS(J,I), FTRIPS(J,I):',J,ITRIP1+II
C      PRINT *,TRIPS(J,MOD(ITRIP1+II,20000)),
C     C      FTRIPS(J,MOD(ITRIP1+II,20000))
C
  10      CONTINUE
   5  CONTINUE
C      PRINT *,' FTRIPS BEFORE:',(JJ,FTRIPS(1,JJ),JJ=0,19999,20)
C
      RETURN
      END
C
C
C
      SUBROUTINE GENREP (ISTATS,COFINM,NUMCOF)
C
*************************************************************************
*
*      TITLE:  GENERATE REPORT
*
*      FILE NAME:  GENREP.FOR
*
*      PURPOSE:  RECIEVES THE RUN STATISTICS TO GENERATE A REPORT WITH.
*
*      HISTORY:
*
*        AUTHOR               DATE            CHANGE
*        ------               ----            ------
*        S.J.KEMPLER          12/1/88         ORIGINAL CODE
*                             2/16/89         ADD VERSION NUMBER TO OUTPUT
*                             4/10/89         ADD NO. OF COEFFICIENTS TO
*                                             OUTPUT.
*
*      CALLING SEQUENCE:  SUBROUTINE GENREP (ISTATS,COFINM,NUMCOF)
*
*      ARGUMENT LIST:
*
*       NAME             TYPE   USE    PURPOSE
*       ----             ----   ---    -------
*
*        ISTATS(10)        I*4     I     ISTATS(1)=NUMBER OF INTERPOLATED
*                                        TRIPLETS ADDED; ISTATS(2)=NUMBER OF
*                                        INTERPOLATED PERIODS; ISTATS(3)=N0.OF
*                                        FILL TRIPLETS ADDED; ISTATS(4)= NUMBER
*                                        OF GAPS AND LARGE MISSING DATA PERIODS.
*                                        ISTATS(5)=TOTAL NUMBER OF TRIPS.
*	COFINM             C*8     I     CONTAINS ALPHANUMERIC NAME OF
*                                        FILTER COEFFICIENTS FILE USED.
*       NUMCOF             I*4     I     NUMBER OF FILTER COEFICIENTS.
*
*      MODULES REFERENCED:  NONE
*
*      NOTES:  RECEIVES THE STATISTICAL INFORMATION FROM THE RUN AND
*              PRINT ON UNIT 6 (LINE PRINTER).
*
*      PDL:
*
*        CALCULATE PERCENTAGE OF INTERPOLATED TRIPLETS ADDED
*        CALCULATE PERCENTAGE OF NUMBER OF FILL TRIPLETS ADDED
*        PRINT COEFFICIENT FILE NAME USED
*        PRINT INTERPOLATED TRIPLETS INFORMATION
*        PRINT FILL TRIPLETS INFORMATION
*        PRINT NUMBER OF "GAP" PERIODS OF DATA
*        RETURN
*        END
*
*******************************************************************************
C
      INTEGER*4 ISTATS(10)
      CHARACTER*8 COFINM
C
      IF (ISTATS(5).NE.0) THEN
       STINTR = FLOAT(ISTATS(1)) / FLOAT(ISTATS(5)) * 100.0
       STFILL = FLOAT(ISTATS(3)) / FLOAT(ISTATS(5)) * 100.0
      ELSE
       STINTR = 0.0
       STFILL = 0.0
      END IF
C
      WRITE(4,98)
C
  98  FORMAT(' ***************************************************')
C
      WRITE(4,99)
C
  99  FORMAT(' VOYAGER FILTERING PROGRAM, VERSION 2.3, 10 APR 1989')
C
      WRITE(4,100) COFINM
C
 100  FORMAT(' NAME OF COEFFICIENT FILE USED: ',A8)
C
      WRITE(4,106) NUMCOF
C
 106  FORMAT(' NUMBER OF FILTER COEFFICIENTS USED: ',I5)
C
      WRITE(4,103) ISTATS(5)
C
 103  FORMAT(' THE TOTAL NUMBER OF DATA POINTS RECEIVED BY VOYFILT: ',
     F       I6/)
C
      WRITE(4,101) ISTATS(1), STINTR
C
 101  FORMAT(' THE NUMBER OF INTERPOLATED DATA POINTS ADDED:        ',
     F       I6,2X,'(',F6.2,'% )' )
C
      WRITE(4,105) ISTATS(2)
C
 105  FORMAT(' THE NUMBER OF PERIODS OF INTERPOLATED DATA ADDED:    ',
     F       I6)
C
      WRITE(4,102) ISTATS(3), STFILL
C
 102  FORMAT(' THE NUMBER OF FILL DATA POINTS ADDED:                ',
     F       I6,2X,'(',F6.2,'% )' )
C
C  DECREMENT BY ONE DUE TO BEGINNING OF DATA GAP BEING COUNTED
C
      WRITE(4,104) ISTATS(4)-1
 104  FORMAT(' THE NUMBER OF DATA "GAP" PERIODS ENCOUNTERED:        ',
     F       I6)
C
      RETURN
      END
C                                                 [U2RFT.MATH]LSPOL.FOR
C
      SUBROUTINE LSPOL(F, X, NX, N, C, RMS, ERR, A, FF, KKK)
C
C     FIT POLYNOMIAL  P(X)  TO DATA  (X,F)     P(X) approx F
C
C     P(X) = C(1) + C(2)*X + C(3)*X*X + ....
C          = (( .... C(3))*X + C(2))*X + C(1)
C
C                                         (R.F. THOMPSON  MARCH 1976)
C        N  = DEGREE OF POLYNOMIAL
C        NX = NUMBER OF DATA POINTS
C
C    ARRAY   MIN DIMENSION
C        F   NX
C        X   NX
C        C   N+1
C      ERR   NX
C        A   3*NX*(N+1)
C       FF   NX
C      KKK   4*(N+1)
C
C
      DIMENSION F(1), X(1), C(1), ERR(1), A(1), FF(1), KKK(1)
      DIMENSION Y(2)
C
      DATA PIVMI/1.0E-6/
C
C
      N1 = N + 1
      DO 10 I=1,NX
   10 FF(I) = F(I)
C        FORM COEFFICIENT MATRIX  A  FOR SYSTEM OF OVERDETERMINED
C        EQUATIONS   A*C = F  .
      IJ = 0
      DO 4 I=1,NX
      EX = X(I)
      IJ = IJ + 1
      A(IJ) = 1.0
      DO 4 J=2,N1
      IJ = IJ + 1
    4 A(IJ) = A(IJ-1)*EX
C
C
      I2 = 1 +   NX*(N+1)
      I3 = 1 + 2*NX*(N+1)
C
      K2 =   N + 2
      K3 = 2*N + 3
      K4 = 3*N + 4
C
C        COMPUTE LEAST SQUARES  N-TH  DEGREE POLYNOMIAL
      IRANK = 0
      CALL SEUDO(A, C, FF, NX, N1, PIVMI,IRANK,
     C  KKK(1),KKK(K2),KKK(K3),KKK(K4), A(I2), A(I3), Y)
C        COMPUTE RMS
      RMS = 0.0
      DO 6 I=1,NX
      EX = X(I)
      PEX = EX*C(N1) + C(N)
      IF(N - 1) 5,5,44
   44 DO 45 JJ=2,N
      J = N1 - JJ
   45 PEX = EX*PEX + C(J)
C
    5 ERR(I) = F(I) - PEX
      RMS = RMS + ERR(I)**2
    6 CONTINUE
C
      RMS = SQRT(RMS/NX)
C
      RETURN
      END


C                S  E  U  D  O



C                                DRA1:[U2RFT.MATH]SEUDO.FOR


C       THIS PROGRAM CAN BE USED TO COMPUTE THE LEAST SQUARES
C       SOLUTION  X  OF AN OVERDETERMINED SYSTEM OF LINEAR
C       EQUATIONS    A*X = B  .



      SUBROUTINE SEUDO(A,X,B,M,N,PIVMI,IRANK,KP1,LP2,KKK,LLL,S,T,Y)

C**** REAL*8 A, X, B, S, T, SUM, BI, AIJ, AIK, AIL  !  DOUBLE PRECISION
      DIMENSION A(1),X(1),B(1),KP1(1),LP2(1),KKK(1),
     C          LLL(1),S(1),T(1),Y(1)
C
C        MATRIX  A  HAS  M  ROWS AND  N  COLUMNS.
C        GIVEN MATRIX  A  AND VECTOR B,
C        THIS PROGRAM COMPUTES A FACTORISATION OF THE PSEUDO-INVERSE
C        OF  A  AND SOLVES  A*X=B  FOR THE MINIMAL NORM LEAST SQUARES
C        SOLUTION.
C     NOTES
C    1  THE INPUT CONTENTS OF ARRAY  A  ARE REPLACED BY THE
C       FACTORISATION OF THE PSEUDOINVERSE OF  A.
C    2  DATA VALUES IN ARRAY  B  ARE DESTROYED
C    3  THE USER MUST SET  IRANK = 0  BEFORE FIRST CALL TO SEUDO.
C       SEVERAL PROBLEMS INVOLVING THE SAME MATRIX  A  CAN BE
C       EFFICIENTLY SOLVED IF THE USER CHANGES ONLY DATA VECTOR  B
C       BETWEEN CALLS TO SEUDO.
C       IRANK = 0          FIRST CALL.  ARRAY  A  CONTAINS MATRIX  A.
C                          SEUDO REPLACES CONTENTS OF ARRAY  A  WITH
C                          FACTORISATION OF PSEUDOINVERSE OF A , AND
C                          SETS IRANK = THE RANK OF A
C       IRANK = RANK OF A  SUCCEEDING CALLS. ARRAY  A  CONTAINS
C                          FACTORIZATION COMPUTED BY FIRST CALL.  USER
C                          MUST CHANGE NOTHING EXCEPT  B  BETWEEN CALLS.
C    4  MINIMUM ARRAY DIMENSIONS:
C                                 A  M*N
C                                 X  N
C                                 B  M
C                               KP1  N
C                               LP2  N
C                               KKK  N
C                               LLL  N
C                                 S  N*N
C                                 T  N*N
C                                 Y  N*N
C
C    5  SET PIVMI = 1.0 E-8  FOR MOST PROBLEMS ON VAX 11/780
C    6  SUBROUTINES NEEDED.... LUDCO, LINSY
C    7   REFERENCE = PETERS AND WILKINSON, THE LEAST SQUARES PROBLEM AND
C        PSEUDO-INVERSES, COMP. J. AUGUST 1970, 309-316 .
C
C                            PROGRAMMED BY R.F. THOMPSON, 4 JULY 1973
C
C
C        MODIFICATIONS
C  29 AUG 74  STATEMENTS 101,102,103,104 PERMIT SECOND AND SUCCESSIVE
C        PROBLEMS WITH SAME TALL MATRIX  A  OF MAXIMUM RANK TO BE
C        SOLVED EFFICIENTLY.
C   7 MAR 75  STATEMENTS  50  THROUGH  60  COMPUTE  S  WITHOUT FIRST
C        STORING MATRIX  A  IN ARRAY  T .  THIS LOWERS MINIMUM DIMENSION
C        OF ARRAY  T  TO  IRANK  WHEN  A  IS TALL OR SQUARE OF MAXIMUM
C        RANK, AND TO  IRANK*IRANK  OTHERWISE.
C   3 APR 75  DOUBLE PRECISION CODE FOR ACCUMULATING L(TRANSPOSE)*B
C        AND L(TRANSPOSE)*L SIGNIFICANTLY REDUCES ERROR BOUND OF
C        COMPUTED SOLUTION.
C
C
      JUMP=2
      IF(IRANK) 101,101,102
  101 JUMP=1
C
C        FIRST, USE COMPLETE PIVOTING TO COMPUTE  RANK(A)=IRANK, AND
C        FACTOR  P1*A*P2=LT*UT, WHERE P1,P2 ARE PERMUTATION MATRICES
C        GIVING THE ROW AND COLUMN INTERCHANGES RESPECTIVELY.
C        LT  IS  M BY IRANK , UNIT LOWER TRAPEZOIDAL, OF RANK IRANK
C        UT  IS  IRANK BY N , UPPER TRAPEZOIDAL, OF RANK IRANK
C        THE PERMUTATIONS  P1 , P2  ARE TABULATED IN  KP1 , LP2  RESP.
C        THE ALGORITHM REPLACES  A  WITH THE ESSENTIAL PARTS OF LT, UT .
      MR=N
      IF(M-N) 1,2,2
    1 MR=M
C        MR=MAXIMUM RANK
C        NS=STEP NUMBER
    2 DO 14 NS=1,MR
C        DETERMINE PIVOT, A(IJ) OF MAXIMUM MODULUS, (I,J) .GE. (NS,NS)
      PIV=PIVMI
      DO 4 I=NS,M
      DO 4 J=NS,N
      IJ=J + N*(I-1)
      IF(PIV - ABS(A(IJ))) 3,4,4
    3 APIV=A(IJ)
      PIV=ABS(APIV)
      KP1(NS)=I
      LP2(NS)=J
    4 CONTINUE
      IF(PIV - PIVMI) 18,18,5
C        INTERCHANGE ROWS  NS  AND  KP1(NS)
    5 IF(NS-KP1(NS)) 6,8,6
    6 N1=N*(NS-1)
      K1=N*(KP1(NS)-1)
      DO 7 J=1,N
      NJ=N1 + J
      KJ=K1 + J
      SAVE=A(NJ)
      A(NJ)=A(KJ)
    7 A(KJ)=SAVE
C        INTERCHANGE COLUMNS  NS  AND  LP2(NS)
    8 IF(NS-LP2(NS)) 9,11,9
    9 DO 10 I=1,M
      NI=N*(I-1)
      IN=NS + NI
      IL=LP2(NS) + NI
      SAVE=A(IN)
      A(IN)=A(IL)
   10 A(IL)=SAVE
   11 IF(NS-MR) 12,15,15
C        COMPUTE ELEMENTS  M(IS) OF  LT  AND  A(IJ)  OF  UT
   12 NS1=NS+1
      DO 13 I=NS1,M
      N1=N*(I-1)
      IS=NS+N1
      A(IS)=A(IS)/APIV
      DO 13 J=NS1,N
      IJ=J + N1
      NSJ=J + N*(NS-1)
   13 A(IJ)=A(IJ) - A(IS)*A(NSJ)
   14 CONTINUE
C
C        MAXIMUM RANK   (BRANCH FROM 11)
   15 IRANK=MR
      IF(N-M) 16,19,19
C        MATRIX  A  IS TALL AND HAS FULL RANK, MR.
C        COMPUTE COLUMN  MR  OF  LT .
   16 MR1=MR+1
      DO 17 I=MR1,M
      IS=MR+N*(I-1)
   17 A(IS)=A(IS)/APIV
      GO TO 19
   18 IRANK=NS-1
   19 CONTINUE
C
C        CODE INVALID FOR IRANK .LE. 1
      IF(IRANK-1) 20,20,22
C  20 WRITE(6,21) IRANK
   20 WRITE(3,21) IRANK
   21 FORMAT(/1X,'***** SEUDO TERMINATED.  RANK ='I3,' .LE. 1')
      RETURN
   22 CONTINUE
C        USE PERMUTATIONS TABULATED IN  KP1  TO PERMUTE  B  TO P1*B.
C        STORE IN B.
  102 DO 23 NS=1,IRANK
      KNS=KP1(NS)
      SAVE=B(NS)
      B(NS)=B(KNS)
   23 B(KNS)=SAVE
C        COMPUTE VECTOR  W = (L TRANSPOSE)*P1*B  (IRANK BY 1)
C        (INNER PRODUCTS OF COLUMNS 1 THROUGH IRANK OF A, WITH B)
C        STORE IN B
      DO 26 J=1,IRANK
      IF(J-M) 24,26,26
   24 I1=J+1
      SUM = B(J)
      DO 25 I=I1,M
      IJ=J+N*(I-1)
      AIJ = A(IJ)
      BI  = B(I)
   25 SUM = SUM + AIJ*BI
      B(J) = SUM
   26 CONTINUE
      GO TO(105,103),JUMP
C        WHEN SOLVING SECOND AND SUCCEEDING PROBLEMS FOR MATRIX  A ,
C        (JUMP=2), TEST IF MATRIX  A  IS TALL AND HAS MAXIMUM RANK,
C        AND EXECUTE APPROPRIATE CODE.
  103 IF(N - M) 104,104,106
  104 IF(IRANK - MR) 106,107,106
C
C
  105 CONTINUE
C        COMPUTE S = (LT TRANSPOSE)*LT      (IRANK BY IRANK)
C        (INNER PRODUCTS OF COLUMNS OF LT)
C             7 MAR 75  NEW CODE, STATEMENTS 50 THROUGH 60,
C             COMPUTE  S  WITHOUT FIRST CREATING EXPLICIT IMAGE
C             OF  LT  IN ARRAY  T  .  THIS REDUCES MINIMUM SIZE OF
C             T  FROM  M BY IRANK  TO  IRANK (WHEN  A  HAS MAXIMUM RANK
C             AND IS TALL OR SQUARE), OR TO  IRANK**2  (OTHERWISE).
      L = IRANK + 1
   50 L = L - 1
      K = L
      GO TO 52
   51 K = K - 1
   52 CONTINUE
C
C        COMPUTE  S(LK) = S(KL) = INNER PRODUCT OF COLUMNS  K  AND  L
C        OF MATRIX  LT.  LT  IS  M BY IRANK  AND IS ZERO ABOVE THE
C        DIAGONAL,  ONE  ON THE DIAGONAL, AND ELEMENTS BELOW
C        THE DIAGONAL ARE STORED IN MATRIX  A.
      IF(K-L) 53,54,54
   53 LK = K + (L-1)*N
      SUM = A(LK)
      GO TO 55
   54 SUM = 1.0
   55 IF(L - M) 56,58,58
   56 I = L
   57 I = I + 1
      II = (I-1)*N
      IK = II + K
      IL = II + L
      AIK = A(IK)
      AIL = A(IL)
      SUM = SUM + AIK*AIL
      IF(I - M) 57,58,58
   58 KL = L + (K-1)*IRANK
      S(KL) = SUM
      LK = K + (L-1)*IRANK
      S(LK) = SUM
      IF(L-1) 60,60,59
   59 IF(K-1) 50,50,51
   60 CONTINUE
C
C
C        GO TO  43  IF MATRIX  A  IS TALL AND HAS MAXIMUM RANK.
      IF(N-M) 30,30,31
   30 IF(IRANK - MR) 31,43,31
   31 CONTINUE
C        COMPUTE T=UT*(UT TRANSPOSE)  (IRANK BY IRANK)
C        (INNER PRODUCTS OF ROWS OF A)
      DO 37 I=1,IRANK
      DO 37 J=I,IRANK
      IJ=J+IRANK*(I-1)
      JI=I+IRANK*(J-1)
      T(IJ)=0.0
      DO 36 KK=J,N
      IK=KK+N*(I-1)
      JK=KK+N*(J-1)
   36 T(IJ)=T(IJ) + A(IK)*A(JK)
   37 T(JI)=T(IJ)
C        FORM  Y = (LT TRANSPOSE)*LT*UT*(UT TRANSPOSE) = S*T
      DO 38 I=1,IRANK
      DO 38 J=1,IRANK
      IJ=J+IRANK*(I-1)
      Y(IJ)=0.0
      DO 38 KK=1,IRANK
      IK=KK+IRANK*(I-1)
      KJ=J+IRANK*(KK-1)
   38 Y(IJ)=Y(IJ) + S(IK)*T(KJ)
C        SOLVE Y*V=W FOR V. STORE IN S
      IP=-1
      PIVM=1.0E-8
      CALL LUDCO(Y,KKK,LLL,IRANK,PIVM,IP,AN,RH)
  106 CALL LINSY(Y,S,B,KKK,LLL,IRANK)
C        COMPUTE Z=(UT TRANSPOSE)*V   (STORE IN X)
C        (INNER PRODUCTS OF COLUMNS OF UT WITH S)
      DO 41 J=1,N
      X(J)=0.0
      ILAST=J
      IF(J-IRANK) 40,40,39
   39 ILAST=IRANK
   40 DO 41 I=1,ILAST
      IJ=J+N*(I-1)
   41 X(J)=X(J) + A(IJ)*S(I)
      GO TO 46
C
C
   43 CONTINUE
C **********TEST************TEST************TEST*************TEST
C     WRITE(6,1201)
C**** WRITE(3,1201)
 1201 FORMAT(/10X,'CODE EXECUTED FOR TALL MATRIX OF MAXIMUM RANK')
C ***************************************************************
C
C        SPECIAL CODE FOR TALL (OR SQUARE) MATRIX OF MAXIMUM RANK, MR.  1 MAY 74
C        AVOIDS SQUARING CONDITION NUMBER.  SEE REFERENCE, FORMULA(8.1)
C        SOLVE (LT TRANSPOSE)*LT*T = W  FOR  T .
      IP = -1
      PIVM = 1.0E-8
      CALL LUDCO(S, KKK, LLL, IRANK, PIVM, IP, AN, RH)
  107 CALL LINSY(S, T, B, KKK, LLL, IRANK)
C
C        SOLVE UT*Z = T  FOR  Z  (BACK SUBSTITUTION)
C        STORE IN  X .
      IJ = MR + N*(MR-1)
      X(MR) = T(MR)/A(IJ)
C        BACK SUBSTITUTION
      MR1 = MR -1
      DO 45 K=1,MR1
      I = MR - K
      SUM = 0.0
      NI1 = N*(I-1)
      J1 = I+1
      DO 44 J=J1,MR
      IJ = J + NI1
   44 SUM = SUM + A(IJ)*X(J)
      II = I + NI1
   45 X(I) = (T(I) - SUM)/A(II)
C        COMPUTE SOLUTION X=P2*Z
C        (USE PERMUTATIONS TABULATED IN LP2 TO PERMUTE ROWS OF X)
C     (N.B.  PREMULTIPLICATION BY  P2  REVERSES ORDER OF PERMUTATIONS)
   46 DO 42 I=1,IRANK
      J=IRANK - I + 1
      LJ=LP2(J)
      SAVE=X(LJ)
      X(LJ)=X(J)
   42 X(J)=SAVE
      RETURN
      END
      SUBROUTINE LINSY (A,X,B,K,L,N)
C
C          THIS SUBROUTINE SOLVES THE LINEAR EQUATION AX=B.
C          A IS AN N BY N MATRIX IN LU FORM (SEE SUBROUTINE LUDCO).
C          THE FIRST N ELEMENTS OF A CONTAIN THE FIRST ROW OF THE MATRIX
C          THE SECOND N ELEMENTS OF ARRAY A CONTAIN THE SECOND ROW OF
C          THE MATRIX, ETC.
C          DURING COMPUTATION, THE CONTENTS OF ARRAY B ARE ESSENTIALLY
C          REPLACED BY (INVERSE OF L)*B, WHICH IS REFERRED TO AS AUGMENT
C          USING THIS AUGMENTED FORM OF B, X IS COMPUTED AS THE
C          SOLUTION TO AX=B BY SIMPLE BACK SUBSTITUTION.
C          AFTER EXECUTION, X WILL CONTAIN THE DIRECT SOLUTION OF
C          AX=B.
C**** REAL*8 A, X, B  !  DOUBLE PRECISION
      DIMENSION A(1),X(1),B(1)
      DIMENSION K(1), L(1)
      N1=N-1
C
C          AUGMENTATION
C
      DO 16 I=1,N1
      MM=L(I)+N*(K(I)-1)
      LL=L(I)
      KK=K(I)
      B(KK)=-B(KK)/A(MM)
      DO 15 J=1,N
      DO 14 JJ=1,I
      IF(K(JJ)-J) 14,15,14
   14 CONTINUE
      M=LL+N*(J-1)
      B(J)= B(J)+B(KK)*A(M)
   15 CONTINUE
   16 CONTINUE
C
C          BACK SUBSTITUTION
C
      MM=L(N)+N*(K(N)-1)
      KK=K(N)
      LL = L(N)
      X(LL)= B(KK)/A(MM)
      DO 18 I=1,N1
      II=N-I
      KK=K(II)
      LL=L(II)
      X(LL)=-B(KK)
      JJ=I
      DO 17 J=1,I
      JJ=JJ-1
      MK=N*(KK-1)
      NN=N-JJ
      LLL=L(NN)
      MM=MK+LLL
   17 X(LL)= X(LL)+A(MM)*X(LLL)
   18 CONTINUE
      RETURN
      END
      SUBROUTINE
     C     LUDCO (A,K,L,N,PIVMI ,IFLAG,ANORM,RHO)
C
C
C
C        THIS SUBROUTINE EFFECTIVELY COMPUTES THE LU DECOMPOSITION
C        OF MATRIX A BY GAUSSIAN ELIMINATION IN AN EXCHANGE ALGORITHM
C        WHICH EMPLOYS EITHER PARTIAL OR COMPLETE PIVOTING.
C        THIS SUBROUTINE COMPUTES THE FOLLOWING
C          1. THE INFINITY NORM, ANORM, OF MATRIX A. (REF. 1 PAGE 98)
C          2. RHO, A PARAMETER WHICH CAN BE USED TO ESTIMATE THE ROUND O
C             ERRORS AND FLOATING POINT ARITHMETIC ERRORS INTRODUCED
C             BY GAUSSIAN ELIMINATION (REF.1 PAGE 102)
C        WHEN CALLING THIS SUBROUTINE....
C          1. A(1) TO A(N) CONTAIN THE FIRST ROW OF MATRIX A,
C             A(N+1) TO A(2N) CONTAIN THE SECOND ROW OF MATRIX A, ETC.
C             HENCE ARRAY A HAS AT LEAST DIMENSION N**2.
C          2. ARRAYS K AND L WILL CONTAIN PIVOTING INFORMATION
C             K(I) CONTAINS THE NUMBER OF THE ITH PIVOT ROW
C             L(I) CONTAINS THE NUMBER OF THE ITH PIVOT COLUMN.
C             THE CALLING PROGRAM MUST DIMENSION K AND L AT LEAST N.
C        UPON RETURN TO THE CALLING PROGRAM, THE LU DECOMPOSED MATRIX A
C        WILL RESIDE IN ARRAY A AND WILL BE IN THE PROPER FORM TO COMPUT
C        THE SOLUTION OF AX=B BY CALLING LINSY.
C        REFERENCES
C          1. FORSYTHE AND MOLER, COMPUTER SOLUTION OF LINEAR ALGEBRAIC
C             SYSTEMS
C          2. STIEFEL, INTRODUCTION TO NUMERICAL ANALYSIS
C
C**** REAL*8 A  !  DOUBLE PRECISION
      DIMENSION A(1),K(1),L(1)
C
C
      LEAP=2
      IF(IFLAG+2) 36,35,36
   35 LEAP=1
   36 CONTINUE
      DO1 I=1,N
      K(I)=0
    1 L(I)=0
      RHO=0.0
      NPIV=0
      GO TO 20
C
C          SEARCH FOR MAXIMUM PIVOT
C
    2 NPIV=NPIV+1
      PIV=0.0
      GO TO (30,40),LEAP
C
C          PARTIAL PIVOTING
C
   30 LT=NPIV
      DO 32 KK=1,N
      DO 33 I=1,NPIV
      IF(K(I)-KK) 33,32,33
   33 CONTINUE
      M=LT+(KK-1)*N
      IF(PIV-ABS(A(M))) 31,32,32
   31 PIV=ABS(A(M))
      KT=KK
   32 CONTINUE
      GO TO 42
C
C          COMPLETE PIVOTING
C
   40 DO 6 KK=1,N
      DO 3 I=1,NPIV
      IF(K(I)-KK) 3,6,3
    3 CONTINUE
      MM=(KK-1)*N
      DO 5 LL=1,N
      DO4 I=1,NPIV
      IF(L(I)-LL) 4,5,4
    4 CONTINUE
      M=MM+LL
      IF(PIV-ABS(A(M))) 41,5,5
   41 PIV=ABS(A(M))
      KT=KK
      LT=LL
    5 CONTINUE
    6 CONTINUE
   42 K(NPIV)=KT
      L(NPIV)=LT
      IF(PIV-RHO) 62,62,61
   61 RHO=PIV
   62 IF(PIV-PIVMI) 63,63,7
   63 CONTINUE
C
C          EXIT, SYSTEM ILL CONDITIONED
C
      IFLAG=NPIV-1
      GO TO (37,38),LEAP
   37 WRITE(3,18) IFLAG
   18 FORMAT(5X,'**** PARTIAL PIVOTING.  ONLY'I5,' EXCHANGES PERFORMED')
      RETURN
   38 WRITE(3,19) IFLAG
   19 FORMAT(5X,'**** THE RANK OF A IS APPROXIMATELY' I5)
      RETURN
C
C
    7 CONTINUE
      IF(NPIV-N) 72,71,72
   71 RHO=RHO/ANORM
      RETURN
C
C          EXCHANGE OF REDUCED MATRIX
   72 CONTINUE
      MK=N*(K(NPIV)-1)
      MP=MK+L(NPIV)
      DO 9 LL=1,N
      DO 8 I=1,NPIV
      IF(LL-L(I)) 8,9,8
    8 CONTINUE
      M=MK+LL
      A(M)=-A(M)/A(MP)
    9 CONTINUE
      DO 12 KK=1,N
      DO 10 I=1,NPIV
      IF(KK-K(I)) 10,12,10
   10 CONTINUE
      M=N*(KK-1)
      KPC=M+L(NPIV)
      APC=A(KPC)
      DO 11 LL=1,N
      DO111 I=1,NPIV
      IF(L(I)-LL) 111,11,111
  111 CONTINUE
      NM=M+LL
      MKLL=MK+LL
      A(NM)=A(NM)+APC*A(MKLL)
   11 CONTINUE
   12 CONTINUE
      GO TO 2
C
C          COMPUTE INFINITY NORM OF A
   20      CONTINUE
      ANORM=0.0
      SUM=0.0
      DO 22 NROW=1,N
      MK=N*(NROW-1)
      DO 21 NCOL=1,N
      M=MK+NCOL
   21 SUM=SUM+ABS(A(M))
      IF(SUM-ANORM) 212,212,211
  211 ANORM=SUM
  212 CONTINUE
      SUM=0.0
   22 CONTINUE
      GO TO 2
      END
