      SUBROUTINE SECGPS (IUNIT,IUO,IOBS,B,NX,G,FATAL)


*** FORM OBS. EQ FROM G FORMAT

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL FULL,FATAL,GETPRM,GETIVF
      LOGICAL LBB,LGF,LCS,LVD,LVA,LVZ,LVS,
     &        LVR,LVG,LVC,LIS,LPS,LPG,LDR,LOS,LAP
      LOGICAL LEB,LLB,LEG,LLG
      CHARACTER*80 CARD
      CHARACTER*1 ID,TC
      DIMENSION B(*),NX(*),G(*)
      DIMENSION ICM((NVECS+1)*15+1)
      DIMENSION KINDS(3*NVECS),ISNS(3*NVECS),JSNS(3*NVECS),LOBS(3*NVECS)
      COMMON /OPRINT/ CRIT,LBB,LGF,LCS,LVD,LVA,LVZ,LVS,
     &                LVR,LVG,LVC,LIS,LPS,LPG,LDR,LOS,LAP
      COMMON /GPS/ MAXVEC
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD
      COMMON /PECHO/ LEB,LLB,LEG,LLG
 
*** INITIALIZE GPS VECTOR TRANSFORMATIONS TO NAD83

      CALL INIT83
 
*** TEST FOR MAXIMUM NUMBER OF VECTORS ALLOWED IN GROUPS

      IF (MAXVEC.GT.NVECS) THEN
        WRITE (6,222) MAXVEC,NVECS
 222    FORMAT ('0ERROR - THE NUMBER OF VECTORS IN A GROUP = ',I2,
     *          ' HAS EXCEEDED THE MAXIMUM ALLOCATED = ',I2,/,
     *  ' *** FATAL -- INCREASE NVECS IN PARAMETER STATEMENTS ***'/)
        CALL ABORT2
      ENDIF

      IF (LGF) THEN
        CALL HEAD
        CALL LINE (2)
        WRITE (6,3)
    3   FORMAT (' ******** GPS OBSERVATIONS *************',/)
      ENDIF

      FULL = .FALSE.
      CALL NEWICM (NICM)

*** ENTER PROCESSING LOOP--EXIT ON END OF FILE

  100 READ (IUNIT,1,END = 777) CARD
    1 FORMAT (A80)
      ID = CARD(1:1)
      IF (.NOT.LEG) THEN
        IF (.NOT.LLG) THEN
          IF (LGF) THEN
            CALL LINE (1)
            WRITE (6,2) CARD
    2       FORMAT (10X,A80)
          ENDIF
        ENDIF
      ENDIF

*** GROUP HEADER RECORD

      IF (ID.EQ.'B') THEN
        IF (FULL) CALL OBSEQW (IUO,FULL,G,B,NX,NVEC,NR,NC,LENG,NICM,ICM,
     &              KINDS,ISNS,JSNS,LOBS,IAUX,IVF,FATAL,IOBS,ITIME)
        READ (CARD,10) IYR1,IMO1,IDY1,IHR1,IMN1,
     &                 IYR2,IMO2,IDY2,IHR2,IMN2,NVEC,I83
   10   FORMAT (1X, 2(I4,4I2),I2,24X,I2)
        TC = 'Z'
        ICODE = 25

        IF (NVEC.LT.1 .OR. NVEC.GT.MAXVEC) THEN
          WRITE (6,11) CARD,NVEC,MAXVEC
   11     FORMAT (/,' ****',/,A80,
     &            /,5X,'ERROR - NVEC = ',I5,' EXCEEDS MAXVEC=',I5)
          CALL ABORT2
        ENDIF
        CALL TOMNT (IYR1,IMO1,IDY1,IHR1,IMN1,TC,IOLD)
        CALL TOMNT (IYR2,IMO2,IDY2,IHR2,IMN2,TC,INEW)
        ITIME = (IOLD + INEW) / 2
        IF ( .NOT.GETPRM(ICODE,ITIME,IAUX) ) IAUX = 0
        IF ( .NOT.GETIVF(ICODE,ITIME,IVF) ) IVF = 0
        IF (I83.GE.7) I83=5
        NR = 3 * NVEC
        LENG = (NVEC + 1) * 3 + 1
        IF (NCD.GT.0) LENG = LENG + 12 * (NVEC + 1)
        NC = NR + 5 + LENG
        CALL NEWGRP (NVEC,IUNIT,G,NR,NC,FULL,ICM,NICM,KINDS,
     &               ISNS,JSNS,LOBS,IOBS,IAUX,B,ITIME,I83)

*** EXTRA MEMBER RECORDS

      ELSEIF (ID.EQ.'C') THEN
        WRITE (6,20) CARD
   20   FORMAT ('0ERROR - TOO MANY C RECORDS',/,
     &          ' BAD G-FILE STRUCTURE--',A80)
        CALL ABORT2

*** CORRELATION RECORDS

      ELSEIF (ID.EQ.'D') THEN
        CALL LOADCR (G,CARD,NR,NC)
      ENDIF

      GO TO 100

*** END OF PROCESSING--END OF FILE ENCOUNTERED

  777 IF (FULL) CALL OBSEQW (IUO,FULL,G,B,NX,NVEC,NR,NC,LENG,NICM,ICM,
     &              KINDS,ISNS,JSNS,LOBS,IAUX,IVF,FATAL,IOBS,ITIME)
      IF (LGF) THEN
        CALL LINE (2)
        WRITE (6,4)
    4   FORMAT ('0*********** END OF GPS OBSERVATIONS *********')
      ENDIF

      RETURN
      END
      SUBROUTINE INIT83

*** INITIALIZE GPS VECTOR TRANSFORMATIONS TO NAD83

      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      IMPLICIT INTEGER (I-N)

      COMMON /CONST/  PI, PI2, RAD
      COMMON /XLATE/  RX(6), RY(6), RZ(6), SCALE(6),
     +                DELX(6), DELY(6), DELZ(6), COSRX(6), SINRX(6),
     +                COSRY(6), SINRY(6), COSRZ(6), SINRZ(6)

*** GPS GFILE CODE '01' WGS 72 TO NAD83

      RX(1) = 0.0D0/3600.D0/RAD
      COSRX(1) = DCOS(RX(1))
      SINRX(1) = DSIN(RX(1))
      RY(1) = 0.0D0/3600.D0/RAD
      COSRY(1) = DCOS(RY(1))
      SINRY(1) = DSIN(RY(1))
      RZ(1) = -0.554D0/3600.D0/RAD
      COSRZ(1) = DCOS(RZ(1))
      SINRZ(1) = DSIN(RZ(1))
      SCALE(1) = +0.2263D0/1.0D6
      DELX(1) = 0.0D0
      DELY(1) = 0.0D0
      DELZ(1) = +4.5D0

*** GPS GFILE CODE '02' WGS 84 TO NAD83

      RX(2) = 0.D0
      COSRX(2) = 1.D0
      SINRX(2) = 0.D0
      RY(2) = 0.D0
      COSRY(2) = 1.D0
      SINRY(2) = 0.D0
      RZ(2) = 0.D0
      COSRZ(2) = 1.D0
      SINRZ(2) = 0.D0
      SCALE(2) = 0.D0
      DELX(2) = 0.D0
      DELY(2) = 0.D0
      DELZ(2) = 0.D0

*** GPS GFILE CODE '03' WGS 72 TO NAD83
 
      RX(3) = 0.0D0/3600.D0/RAD
      COSRX(3) = DCOS(RX(3))
      SINRX(3) = DSIN(RX(3))
      RY(3) = 0.0D0/3600.D0/RAD
      COSRY(3) = DCOS(RY(3))
      SINRY(3) = DSIN(RY(3))
      RZ(3) = -0.554D0/3600.D0/RAD
      COSRZ(3) = DCOS(RZ(3))
      SINRZ(3) = DSIN(RZ(3))
      SCALE(3) = +0.2263D0/1.0D6
      DELX(3) = 0.0D0
      DELY(3) = 0.0D0
      DELZ(3) = +4.5D0
 
*** GPS GFILE CODE '04' WGS 84 TO NAD83
 
      RX(4) = 0.D0
      COSRX(4) = 1.D0
      SINRX(4) = 0.D0
      RY(4) = 0.D0
      COSRY(4) = 1.D0
      SINRY(4) = 0.D0
      RZ(4) = 0.D0
      COSRZ(4) = 1.D0
      SINRZ(4) = 0.D0
      SCALE(4) = 0.D0
      DELX(4) = 0.D0
      DELY(4) = 0.D0
      DELZ(4) = 0.D0
 
*** GPS GFILE CODE '05' ITRF(VLBI)1989 TO NAD83
 
      RX(5) = 0.0275D0/3600.D0/RAD
      COSRX(5) = DCOS(RX(5))
      SINRX(5) = DSIN(RX(5))
      RY(5) = 0.0155D0/3600.D0/RAD
      COSRY(5) = DCOS(RY(5))
      SINRY(5) = DSIN(RY(5))
      RZ(5) = 0.0107D0/3600.D0/RAD
      COSRZ(5) = DCOS(RZ(5))
      SINRZ(5) = DSIN(RZ(5))
      SCALE(5) = 0.0D0/1.0D6
      DELX(5) = 0.9191D0
      DELY(5) = -2.0182D0
      DELZ(5) = -0.4835D0
 
*** GPS GFILE CODE '06' NEOS 1990 (PROVISIONAL) TO NAD83

       RX(6) = 0.0229D0/3600.D0/RAD
      COSRX(6) = DCOS(RX(6))
      SINRX(6) = DSIN(RX(6))
      RY(6) = 0.0249D0/3600.D0/RAD
      COSRY(6) = DCOS(RY(6))
      SINRY(6) = DSIN(RY(6))
      RZ(6) = 0.0099D0/3600.D0/RAD
      COSRZ(6) = DCOS(RZ(6))
      SINRZ(6) = DSIN(RZ(6))
      SCALE(6) = 0.0D0/1.0D6
      DELX(6) = 0.8845D0
      DELY(6) = -2.0399D0
      DELZ(6) = -0.4835D0
 
      RETURN
      END
      SUBROUTINE NEWGRP (NVEC,IUNIT,G,NR,NC,FULL,ICM,NICM,KINDS,
     &                   ISNS,JSNS,LOBS,IOBS,IAUX,B,ITIME,I83)

*** LOAD AN OBSERVATION GROUP INTO WORK SPACE

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL FULL,LBV,GETSSN
      LOGICAL LEB,LLB,LEG,LLG
      LOGICAL LBB,LGF,LCS,LVD,LVA,LVZ,LVS,
     &        LVR,LVG,LVC,LIS,LPS,LPG,LDR,LOS,LAP
      LOGICAL LMSL,LSS,LUP
      CHARACTER*80 CARD
      CHARACTER*1 ID,CODE
      DIMENSION G(NR,NC)
      DIMENSION ICM((NVECS+1)*15+1)
      DIMENSION KINDS(3*NVECS),ISNS(3*NVECS),JSNS(3*NVECS),LOBS(3*NVECS)
      DIMENSION IC(31),B(*)
      COMMON /PECHO/ LEB,LLB,LEG,LLG
      COMMON /ECHO/ VSD,LBV
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD
      COMMON /OPRINT/ CRIT,LBB,LGF,LCS,LVD,LVA,LVZ,LVS,
     &                LVR,LVG,LVC,LIS,LPS,LPG,LDR,LOS,LAP
      COMMON /STATCT/ N84,N85,N86,N89,NDIR,NANG,NGPS,NZD,NDS,NAZ,NQQ,
     &                NREJ,NGPSR
      COMMON /OPT/ AX,E2,DMSL,DGH,VM,VP,CTOL,ITMAX,ITMIN,IMODE,
     &             LMSL,LSS,LUP
      COMMON /VF/ VFACTR(31)

*** CLEAR WORK AREA

      DO 1 I = 1,NR
        DO 2 J = 1,NR
          G(I,J) = 0.D0
    2   CONTINUE
    1 CONTINUE

*** LOOP OVER NVEC MEMBER RECORDS  ('C')

      DO 10 I = 1,NVEC
        J1 = (I - 1) * 3 + 1
        J2 = J1 + 1
        J3 = J2 + 1
  100   READ (IUNIT,3,END = 666) CARD
    3   FORMAT (A80)
        ID = CARD(1:1)
        IF ( ID.EQ.'A' .OR. ID.EQ.'B' .OR. ID.EQ.'D') GO TO 666
        IF ( ID.EQ.'C' ) THEN
          READ (CARD,4) ISSN,JSSN,DX,SDX,DY,SDY,DZ,SDZ,CODE
    4     FORMAT (1X,I4,   I4,    3(F11.4,F5.4),A1)
        ELSEIF (ID.EQ.'F')THEN
          READ (CARD,44) ISSN,JSSN,DX,SDX,DY,SDY,DZ,SDZ,CODE
   44     FORMAT( 1X,I4,I4,3(F13.4,f5.4),A1)
        ELSE
          GO TO 100
        ENDIF

          IF ( .NOT.GETSSN(ISSN,ISN) ) THEN
            CALL LINE (3)
            WRITE (6,7) CARD
   7        FORMAT ('0ERROR - NO *80* RECORD FOR--',A80,/)
            CALL ABORT2
          ELSEIF ( .NOT.GETSSN(JSSN,JSN) ) THEN
            CALL LINE (3)
            WRITE (6,7) CARD
            CALL ABORT2
          ENDIF

*** CONVERT GPS VECTOR TO DATUM (NAD83)

          IF (I83 .NE. 0) CALL TO83 (DX,DY,DZ,ISN,I83,B,ITIME)

          ISNS(J1) = ISN
          ISNS(J2) = ISN
          ISNS(J3) = ISN

          JSNS(J1) = JSN
          JSNS(J2) = JSN
          JSNS(J3) = JSN

          KINDS(J1) = 4
          KINDS(J2) = 5
          KINDS(J3) = 6

          G(J1,NC) = DX
          G(J2,NC) = DY
          G(J3,NC) = DZ

          IF (CODE.EQ.' ') THEN
            G(J1,J1) = SDX * DSQRT(VFACTR(4))
            G(J2,J2) = SDY * DSQRT(VFACTR(5))
            G(J3,J3) = SDZ * DSQRT(VFACTR(6))

            IOBS = IOBS + 1
            LOBS(J1) = IOBS
            IOBS = IOBS + 1
            LOBS(J2) = IOBS
            IOBS = IOBS + 1
            LOBS(J3) = IOBS
            NOBS = NOBS + 3
            NGPS = NGPS + 1

*** IF REJECTED VECTOR, THEN USE WEIGHTS OF 100 METERS

          ELSE
            G(J1,J1) = 100.D0
            G(J2,J2) = 100.D0
            G(J3,J3) = 100.D0
            LOBS(J1) = IOBS
            LOBS(J2) = IOBS
            LOBS(J3) = IOBS
            NREJ = NREJ + 3
            NGPSR = NGPSR + 1
          ENDIF

*** ECHO GPS OBSERVATION

          IF (.NOT.LLG) THEN
            IF (CODE.EQ.' ') THEN
              CALL LINE (1)
              WRITE (6,5) IOBS,CARD
    5         FORMAT (I7,3X,A80)
            ELSE
              CALL LINE (1)
              WRITE (6,11) CARD
   11         FORMAT (10X,A80)
            ENDIF
          ENDIF

          CALL FORMIC (KINDS(J1),ISNS(J1),JSNS(J1),IAUX,IC,L,B)
          CALL PUTICM (IC,L,ICM,NICM)

*** THE FOLLOWING FOUR LINES ARE CURRENTLY REDUNDANT WITH THE PRIOR TWO
*** OF COURSE, IF CHANGES TO THE GPS OBSERVATION MODEL ARE MADE,
*** THEN IT MAY BE NECESSARY TO REMOVE COMMENTS ON THESE FOUR LINES

***       CALL FORMIC (KINDS(J2),ISNS(J2),JSNS(J2),IAUX,IC,L,B)
***       CALL PUTICM (IC,L,ICM,NICM)
***       CALL FORMIC (KINDS(J3),ISNS(J3),JSNS(J3),IAUX,IC,L,B)
***       CALL PUTICM (IC,L,ICM,NICM)

   10 CONTINUE

*** FELL THRU LOOP--NORMAL END OF PROCESSING

      FULL = .TRUE.
      RETURN

*** BAD G-FILE STRUCTURE

  666 WRITE (6,667) CARD
  667 FORMAT ('0ERROR - BAD G-FILE STRUCTURE--',A80)
      CALL ABORT2
      RETURN
      END
      SUBROUTINE TO83 (DX, DY, DZ, ISN, I83, B, ITIME)

*** CONVERT GPS VECTOR TO NAD83

      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION B(*)

*** GET 83 AND CONVERT TO SATELLITE SYSTEM

      CALL GETXYZ(X,Y,Z,B,ISN,ITIME)
      CALL FROM83 (I83, X, Y, Z, XI, YI, ZI)

*** ADD VECTOR IN SATELLITE SYSTEM

      XJ = XI + DX
      YJ = YI + DY
      ZJ = ZI + DZ

*** CONVERT ENDPOINTS BACK TO 83

      CALL GET83 (I83, XI, YI, ZI, X1, Y1, Z1)
      CALL GET83 (I83, XJ, YJ, ZJ, X2, Y2, Z2)

*** GET VECTOR IN 83

      DX = X2 - X1
      DY = Y2 - Y1
      DZ = Z2 - Z1

      RETURN
      END
      SUBROUTINE FROM83 (I83, X83, Y83, Z83, X, Y, Z)

*** CONVERT X,Y,Z FROM ADJUSTMENT DATUM

      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      IMPLICIT INTEGER (I-N)
      COMMON /XLATE/  RX(6), RY(6), RZ(6), SCALE(6),
     +                DELX(6), DELY(6), DELZ(6), COSRX(6), SINRX(6),
     +                COSRY(6), SINRY(6), COSRZ(6), SINRZ(6)

*** DO THE TRANSLATION

      X = X83 - DELX(I83)
      Y = Y83 - DELY(I83)
      Z = Z83 - DELZ(I83)

*** ROTATE ABOUT THE X AXIS

      XRX = X
      YRX = Y*COSRX(I83) - Z*SINRX(I83)
      ZRX = Z*COSRX(I83) + Y*SINRX(I83)

*** ROTATE ABOUT THE Y AXIS
 
      XRXY = XRX*COSRY(I83) + ZRX*SINRY(I83)
      YRXY = YRX
      ZRXY = ZRX*COSRY(I83) - XRX*SINRY(I83)
 
*** ROTATE ABOUT THE Z AXIS
 
      X = XRXY*COSRZ(I83) - YRXY*SINRZ(I83)
      Y = YRXY*COSRZ(I83) + XRXY*SINRZ(I83)
      Z = ZRXY
 
*** APPLY THE SCALE CHANGE
 
      XNEW = X - SCALE(I83)*X
      YNEW = Y - SCALE(I83)*Y
      ZNEW = Z - SCALE(I83)*Z
 
      RETURN
      END
      SUBROUTINE GET83 (I83, XOLD, YOLD, ZOLD, XNEW, YNEW, ZNEW)
        
*** CONVERT X,Y,Z TO ADJUSTMENT DATUM
 
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      IMPLICIT INTEGER (I-N)
 
      COMMON /XLATE/  RX(6), RY(6), RZ(6), SCALE(6),
     +                DELX(6), DELY(6), DELZ(6), COSRX(6), SINRX(6),
     +                COSRY(6), SINRY(6), COSRZ(6), SINRZ(6)
 
*** DO THE TRANSLATION
 
      X = XOLD + DELX(I83)
      Y = YOLD + DELY(I83)
      Z = ZOLD + DELZ(I83)
 
*** ROTATE ABOUT THE X AXIS
 
      XRX = X
      YRX = Y*COSRX(I83) + Z*SINRX(I83)
      ZRX = Z*COSRX(I83) - Y*SINRX(I83)
 
*** ROTATE ABOUT THE Y AXIS
 
      XRXY = XRX*COSRY(I83) - ZRX*SINRY(I83)
      YRXY = YRX
      ZRXY = ZRX*COSRY(I83) + XRX*SINRY(I83)
 
*** ROTATE ABOUT THE Z AXIS
 
      X = XRXY*COSRZ(I83) + YRXY*SINRZ(I83)
      Y = YRXY*COSRZ(I83) - XRXY*SINRZ(I83)
      Z = ZRXY
 
*** APPLY THE SCALE CHANGE
 
      XNEW = X + SCALE(I83)*X
      YNEW = Y + SCALE(I83)*Y
      ZNEW = Z + SCALE(I83)*Z
 
      RETURN
      END

      SUBROUTINE LOADCR (G,CARD,NR,NC)

*** LOAD CORRELATIONS FROM A RECORD

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      CHARACTER*80 CARD
      DIMENSION G(NR,NC)
      DIMENSION IROW(5),ICOL(5),COR(5)

      READ (CARD,1) ( IROW(I),ICOL(I),COR(I), I = 1,5 )
    1 FORMAT (1X, 5(2I3,F9.7) )

      DO 10 I = 1,5
        IR = IROW(I)
        IC = ICOL(I)
        C = COR(I)
        IF (IR.NE.0 .AND. IC.NE.0) THEN
          IF (DABS(C).GT.1.D0) THEN
            WRITE (6,2) IR,IC,C
    2       FORMAT ('0ERROR - ILLEGAL CORRELATION--',2I5,D20.10)
            CALL ABORT2
          ELSEIF (IR.LT.1 .OR. IR.GT.NR .OR.
     &            IC.LT.1 .OR. IC.GT.NR ) THEN
            WRITE (6,3) IR,IC
    3       FORMAT ('0ERROR - BAD CORRELATION INDEX--',2I5)
            CALL ABORT2
          ELSEIF (IR.EQ.IC) THEN
            CALL LINE (1)
            WRITE (6,4) IR,IC,C
    4       FORMAT (' *** WARNING -- DIAGONAL CORRELATION ',2I5,D20.10)
          ELSE
            G(IR,IC) = C
            G(IC,IR) = C
          ENDIF
        ENDIF
   10 CONTINUE

      RETURN
      END
      SUBROUTINE OBSEQW (IUO,FULL,G,B,NX,NVEC,NR,NC,LENG,NICM,ICM,
     &              KINDS,ISNS,JSNS,LOBS,IAUX,IVF,FATAL,IOBS,ITIME)

*** COMPUTE & DECORRELATE OBS.EQ. FROM WORK AREA

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL ADDCON,FULL,FATAL,NOBIGV
      DIMENSION G(NR,NC),B(*),NX(*)
      DIMENSION ICM((NVECS+1)*15+1)
      DIMENSION KINDS(3*NVECS),ISNS(3*NVECS),JSNS(3*NVECS),LOBS(3*NVECS)
      DIMENSION ICX(31),CX(31)
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD

*** ABORT/PRINT LARGE RESIDUALS

      NOBIGV = .FALSE.

*** DEFINE SIGNAL FOR A GPS OBS. EQ. SECTION

      ISIGNL = 1000

*** STORE STANDARD ERROR OF OBSERVATIONS AND
*** CORRELATION IN G-STRUCTURE

      LENG = (NVEC + 1) * 3 + 1
      IF (NCD .GT. 0) LENG = LENG + 12 * (NVEC + 1)
      NC = NR + 5 + LENG
      CALL GLOCAT(N1,N2,N3,N4,N5,NVEC,LENG)
      DO 100 J = 1, NVEC
         I = 3 * (J-1) + 1
         G(I,N5) = G(I,I)
         G(I+1,N5) = G(I+1,I+1)
         G(I+2,N5) = G(I+2,I+2)
         G(I,N5+1) = G(I,I+1)
         G(I+1,N5+1) = G(I+1,I+2)
         G(I+2,N5+1) = G(I,I+2)
  100 CONTINUE 

*** TRANSFORM CORR-TYPE MATRIX TO COVARIANCES

      DO 3 I = 1,NR
        GII = G(I,I)
        DO 8 J = 1,NR
          IF (I.NE.J) G(I,J) = G(I,J) * GII * G(J,J)
    8   CONTINUE
    3 CONTINUE
      DO 5 I = 1,NR
        G(I,I) = G(I,I) * G(I,I)
        G(I,1) = G(I,I)
    5 CONTINUE

*** CHOLESKY FACTOR COVARIANCE MATRIX IN WORK SPACE

      DO 1 I = 1,NR
        I1 = I - 1
        IF (I1.GT.0) THEN
          DO 2 K = 1,I1
            G(I,I) = G(I,I) - G(K,I) * G(K,I)
    2     CONTINUE
        ENDIF

*** SINGULARITY TEST

        IF ( G(I,I).LT.1.D-10 ) THEN
          WRITE (6,7) I
    7     FORMAT ('0ERROR - CORRELATION MATRIX SINGULARITY IN'
     &            ' OBSEQW--',I5)
          CALL ABORT2
        ENDIF

*** RESUME CHOLESKY FACTOR

        G(I,I) = DSQRT( G(I,I) )
        IF ( I + 1 .LE. NR ) THEN
          DO 6 J = I + 1,NR
            IF (I1.GT.0) THEN
              DO 4 K = 1,I1
                G(I,J) = G(I,J) - G(K,I) * G(K,J)
    4         CONTINUE
            ENDIF
            G(I,J) = G(I,J) / G(I,I)
  6       CONTINUE
        ENDIF
    1 CONTINUE

*** UPDATE CONNECTIVITY AND WRITE GPS OBS.EQ. SECTION


      IF ( .NOT. ADDCON(ICM,NICM,NX) ) THEN
        WRITE (6,666) NICM, (ICM(I),I = 1,NICM)
  666   FORMAT ('0ERROR - INSUFFICIENT MEMORY GPS VECTORS',10I3)
        CALL ABORT2
      ENDIF
      WRITE (IUO) ISIGNL,NVEC,IAUX,ICX,CX,LENGX,CMOX,OBX,SDXX,
     &            IOBS,IVF,IAUX,ITIME
      WRITE (IUO) ICM,NICM,KINDS,ISNS,JSNS,LOBS
      CALL GOBSEQ (G,NR,NC,NVEC,LENG,B,ICM,NICM,KINDS,ISNS,JSNS,
     &             LOBS,IAUX,IVF,FATAL,NOBIGV,N4,ITIME,N5)
      DO 10 I = 1,NR
        WRITE (IUO) (G(I,J),J = 1,NC)
   10 CONTINUE

      FULL = .FALSE.
      CALL NEWICM (NICM)

      RETURN
      END
      SUBROUTINE FORMG (IUO,IUO2,NVEC,NR,G,B,FATAL,IVF,IAUX,ITIME)

*** ROUTINE TO RE-COMPUTE OBS.EQ. FOR A GPS CLUSTER

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL FATAL,NOBIGV
      DIMENSION B(*),G(NR,*)
      DIMENSION ICM((NVECS+1)*15+1)
      DIMENSION KINDS(3*NVECS),ISNS(3*NVECS),JSNS(3*NVECS),LOBS(3*NVECS)
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD

*** ABORT/PRINT BIG RESIDUALS

      NOBIGV = .FALSE.

      LENG = (NVEC + 1) * 3 + 1
      IF (NCD.GT.0) LENG = LENG + 12 * (NVEC + 1)
      NC = NR + 5 + LENG

*** READ SUPPORTING INDICIES

      READ (IUO,END=666) ICM,NICM,KINDS,ISNS,JSNS,LOBS
      WRITE (IUO2) ICM,NICM,KINDS,ISNS,JSNS,LOBS

*** LOAD WORK SPACE (G)

      DO 1 I = 1,NR
        READ (IUO,END=666) ( G(I,J),J=1,NC )
    1 CONTINUE

*** COMPUTE AND DECORRELATE OBS. EQ. AND MISCLOSURE

      CALL GOBSEQ (G,NR,NC,NVEC,LENG,B,ICM,NICM,KINDS,ISNS,JSNS,
     &             LOBS,IAUX,IVF,FATAL,NOBIGV,N4,ITIME,N5)

*** UNLOAD WORK SPACE

      DO 500 I = 1,NR
        WRITE (IUO2) ( G(I,J),J = 1,NC )
  500 CONTINUE

      RETURN

*** PREMATURE END OF FILE

  666 WRITE (6,667) NVEC
  667 FORMAT ('0ERROR - PREMATURE FILE END IN FORMG -- NVEC=',I5)
      CALL ABORT2
      RETURN
      END
      SUBROUTINE GOBSEQ (G,NR,NC,NVEC,LENG,B,ICM,NICM,KINDS,ISNS,JSNS,
     &                   LOBS,IAUX,IVF,FATAL,NOBIGV,N4,ITIME,N5)

*** COMPUTE AND DECORRELATE OBS. EQ. AND MISCLOSURE

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL FATAL,NOBIGV
      LOGICAL LMSL,LSS,LUP
      DIMENSION KINDS(3*NVECS),ISNS(3*NVECS),JSNS(3*NVECS),LOBS(3*NVECS)
      DIMENSION ICM((NVECS+1)*15+1)
      DIMENSION B(*),G(NR,NC)
      DIMENSION IC(31),C(31)
      COMMON /OPT/ AX,E2,DMSL,DGH,VM,VP,CTOL,ITMAX,ITMIN,IMODE,
     &             LMSL,LSS,LUP
      COMMON /LASTSN/ ISNX,JSNX,KSNX,ITIMEX,JTIMEX,KTIMEX
      COMMON /LASTXT/ ISNZ,JSNZ,KSNZ,ITIMEZ,JTIMEZ,KTIMEZ
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD

*** COMPUTE THE MAXIUM LENGTH OF COEFF. ARRAYS

      LMAX = 6
      IF (NCD.GT.0) LMAX = LMAX + 24
      IF (IAUX.GT.0) LMAX = LMAX + 1

*** SET LAST STATION ENCOUNTERED TO NULL

      ISNX = 0
      JSNX = 0
      KSNX = 0
      ITIMEX = 0
      JTIMEX = 0
      KTIMEX = 0
      ISNZ = 0
      JSNZ = 0
      KSNZ = 0
      ITIMEZ = 0
      JTIMEZ = 0
      KTIMEZ = 0

*** DETERMINE WORK ARRAY COLUMN POINTERS

      CALL GLOCAT (N1,N2,N3,N4,N5,NVEC,LENG)

*** COMPUTE COEFFICIENTS AND OBSERVATIONS

      DO 10 I = 1,NR
        LENGC = LMAX
        KIND = KINDS(I)
        ISN = ISNS(I)
        JSN = JSNS(I)
        IOBS = LOBS(I)
        CALL COMPOB (KIND,OBS0,B,RDUMMY,ISN,JSN,IAUX,ITIME)
        CALL FORMC (KIND,C,B,ISN,JSN,IAUX,ITIME)
        CALL GETGIC (ISN,JSN,IAUX,ICM,IC,N2,NICM,B)
        CALL COMPIC (IC,C,LENGC)

*** LOAD WORK ARRAY

        CMO = OBS0 - G(I,NC)
        IF (IMODE.EQ.0) CMO = 0.D0
        G(I,N2) = CMO

        DO 2 J = 1,LENG
          K = J + N2
          G(I,K) = 0.D0
    2   CONTINUE

        DO 20 J = 1,LENGC
          G(I,IC(J) ) = C(J)
   20   CONTINUE
   10 CONTINUE

*** DECORRELATE MISCLOSURE

      CALL COMRHS (G,NR)
      CALL PUTRHS (G,NR,NC,N2,N4)

*** TEST FOR LARGE MISCLOSURES

      IF (.NOT.NOBIGV) THEN
        SD = 1.D0
        DO 15 I = 1,NR
          CMO = G(I,N2)
          CALL BIGV (KINDS(I),ISNS(I),JSNS(I),LOBS(I),IVF,CMO,SD,FATAL)
   15   CONTINUE
      ENDIF

*** DECORRELATE COEFFICIENTS

      DO 30 I = 1,LENG
        CALL GETRHS (G,NR,NC,N2,N2+I)
        CALL COMRHS (G,NR)
        CALL PUTRHS (G,NR,NC,N2,N2+I)
   30 CONTINUE

      RETURN
      END
      SUBROUTINE FORMG2 (IUO,IUO2,NVEC,NR,G,B,A,NX,FATAL,
     &                   IVF,IAUX,ITIME)

*** ROUTINE TO RE-COMPUTE OBS.EQ. FOR A GPS CLUSTER

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL FATAL
      DIMENSION B(*),G(NR,*),A(*),NX(*)
      DIMENSION ICM((NVECS+1)*15+1)
      DIMENSION KINDS(3*NVECS),ISNS(3*NVECS),JSNS(3*NVECS),LOBS(3*NVECS)
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD

      LENG = (NVEC + 1) * 3 + 1
      IF (NCD.GT.0) LENG = LENG + 12 * (NVEC + 1)
      NC = NR + 5 + LENG

*** READ SUPPORTING INDICIES

      READ (IUO) ICM,NICM,KINDS,ISNS,JSNS,LOBS
      WRITE (IUO2) ICM,NICM,KINDS,ISNS,JSNS,LOBS

*** LOAD WORK SPACE (G)

      DO 1 I = 1,NR
        READ (IUO) ( G(I,J),J = 1,NC )
    1 CONTINUE

*** COMPUTE AND DECORRELATE OBS. EQ. AND MISCLOSURE

      CALL VFGPS (G,NR,NC,NVEC,LENG,B,ICM,NICM,KINDS,ISNS,JSNS,
     &            LOBS,IAUX,IVF,FATAL,A,NX,N2,N3,N4,ITIME)

*** UNLOAD WORK SPACE

      DO 500 I = 1,NR
        WRITE (IUO2) ( G(I,J),J = 1,NC )
  500 CONTINUE

      RETURN
      END
      SUBROUTINE VFGPS (G,NR,NC,NVEC,LENG,B,ICM,NICM,KINDS,ISNS,JSNS,
     &                  LOBS,IAUX,IVF,FATAL,A,NX,N2,N3,N4,ITIME)

*** COMPUTE AND DECORRELATE OBS. EQ. AND MISCLOSURE

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL FATAL,PROP,FIXVF
      LOGICAL LMSL,LSS,LUP
      DIMENSION KINDS(3*NVECS),ISNS(3*NVECS),JSNS(3*NVECS),LOBS(3*NVECS)
      DIMENSION ICM((NVECS+1)*15+1),CM((NVECS+1)*15+1)
      DIMENSION B(*),G(NR,NC),A(*),NX(*)
      DIMENSION IC(31),C(31)
      COMMON /OPT/ AX,E2,DMSL,DGH,VM,VP,CTOL,ITMAX,ITMIN,IMODE,
     &             LMSL,LSS,LUP
      COMMON /LASTSN/ ISNX,JSNX,KSNX,ITIMEX,JTIMEX,KTIMEX
      COMMON /LASTXT/ ISNZ,JSNZ,KSNZ,ITIMEZ,JTIMEZ,KTIMEZ
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD
      COMMON /VFTB2/ VFS(30),VTV(30),VFRN(30),VFSS(30)

*** LENGTH OF COEFF. ARRAYS

      LENGC = 6
      IF (NCD.GT.0) LENGC = LENGC + 24
      IF (IAUX.GT.0) LENGC = LENGC + 1

*** SET LAST STATION ENCOUNTERED TO NULL

      ISNX = 0
      JSNX = 0
      KSNX = 0
      ITIMEX = 0
      JTIMEX = 0
      KTIMEX = 0
      ISNZ = 0
      JSNZ = 0
      KSNZ = 0
      ITIMEZ = 0
      JTIMEZ = 0
      KTIMEZ = 0

*** DETERMINE WORK ARRAY COLUMN POINTERS

      CALL GLOCAT (N1,N2,N3,N4,N5,NVEC,LENG)

*** (STEP 1.)  COMPUTE COEFFICIENTS

      DO 10 I = 1,NR
        KIND = KINDS(I)
        ISN = ISNS(I)
        JSN = JSNS(I)
        CALL FORMC (KIND,C,B,ISN,JSN,IAUX,ITIME)
        CALL GETGIC (ISN,JSN,IAUX,ICM,IC,N2,NICM,B)

        DO 2 J = 1,LENG
          K = J + N2
          G(I,K) = 0.D0
    2   CONTINUE

        DO 20 J = 1,LENGC
          G(I,IC(J) ) = C(J)
   20   CONTINUE
   10 CONTINUE

*** (STEP 2.)  PROPAGATE FOR STD. DEV. OF RESIDUALS


*** (STEP 3.)  DECORRELATE COEFFICIENTS

      DO 30 I = 1,LENG
        CALL GETRHS (G,NR,NC,N2,N2+I)
        CALL COMRHS (G,NR)
        CALL PUTRHS (G,NR,NC,N2,N2+I)
   30 CONTINUE

*** (STEP 4.)  PROPAGATE FOR REDUNDANCY NUMBER (IF VARIANCE FACTOR)

      IF (IVF.GT.0) THEN
        IF ( .NOT.FIXVF(IVF) ) THEN
          DO 80 IROW = 1,NR

*** BYPASS REJECTION WHEN COMPUTING A REDUNDANCY NUMBER

            IF ( G(IROW,1) .GE. 100.D0 ) GO TO 80
            DO 70 I = 1,NICM
              CM(I) = G(IROW,N2+I)
   70       CONTINUE

            IF ( .NOT.PROP(CM,ICM,NICM,W,A,NX,IFLAG) ) THEN
              IF (IFLAG.EQ.1) THEN
                WRITE (6,668)
  668           FORMAT ('0STATE ERROR IN PROP OF VFGPS')
                CALL ABORT2
              ELSE
                WRITE (6,669)
  669           FORMAT ('0PROFILE ERROR IN PROP OF VFGPS')
                CALL ABORT2
              ENDIF
            ENDIF

            W = W / VFS(IVF)
            RN = 1.D0 - W
            IF ( RN .LT. 1.D-6 ) RN = 0.D0
            VFRN(IVF) = VFRN(IVF) + RN
   80     CONTINUE
        ENDIF
      ENDIF

*** (STEP 5.)  COMPUTE OBSERVATIONS

      DO 90 I = 1,NR
        KIND = KINDS(I)
        ISN = ISNS(I)
        JSN = JSNS(I)
        CALL COMPOB (KIND,OBS0,B,RDUMMY,ISN,JSN,IAUX,ITIME)
        CMO = OBS0 - G(I,NC)
        G(I,N2) = CMO
   90 CONTINUE

*** (STEP 6.)  DECORRELATE MISCLOSURE

      CALL COMRHS (G,NR)
      CALL PUTRHS (G,NR,NC,N2,N4)
      SD = 1.D0
      DO 16 I = 1,NR
        CMO = G(I,N2)
        IF (IVF.GT.0) THEN
          IF ( .NOT.FIXVF(IVF) ) THEN
            VTV(IVF) = VTV(IVF) + CMO * CMO / VFS(IVF)
          ENDIF
        ENDIF
        CALL BIGV (KINDS(I),ISNS(I),JSNS(I),LOBS(I),IVF,CMO,SD,FATAL)
   16 CONTINUE

      RETURN
      END
      SUBROUTINE GLOCAT (N1,N2,N3,N4,N5,NVEC,LENG)

*** DETERMINE WORK ARRAY COLUMN POINTERS

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)

      N1 = 3 * NVEC
      N2 = N1 + 1
      N3 = N2 + LENG
      N4 = N3 + 1
      N5 = N4 + 1

      RETURN
      END
      SUBROUTINE COMRHS (G,N)

*** SINGLE SUBSTITUTION FOR DECORRELATION FROM A CHOL. FACTOR

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION G(N,*)

      N1 = N + 1
      DO 1 I = 1,N
        I1 = I - 1
        TEMP = 0.D0
        IF (I1.GT.0) THEN
          DO 2 K = 1,I1
            TEMP = TEMP + G(K,I) * G(K,N1)
    2     CONTINUE
        ENDIF
        G(I,N1) = (G(I,N1) - TEMP) / G(I,I)
    1 CONTINUE

      RETURN
      END
      SUBROUTINE GETRHS (G,NR,NC,N,M)

*** COPY M-TH COLUMN TO N-TH COLUMN

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION G(NR,NC)

      DO 1 I = 1,NR
        G(I,N) = G(I,M)
    1 CONTINUE

      RETURN
      END
      SUBROUTINE PUTRHS (G,NR,NC,N,M)

*** COPY N-TH COLUMN TO M-TH COLUMN

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION G(NR,NC)

      DO 1 I = 1,NR
        G(I,M) = G(I,N)
    1 CONTINUE

      RETURN
      END
      SUBROUTINE GETGIC (ISN,JSN,IAUX,ICM,IC,N2,NICM,B)

*** GET COLUMN INDEX IN WORK ARRAY FROM STATION NUMBERS

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL GETICM
      DIMENSION ICM((NVECS+1)*15+1)
      DIMENSION IC(31)

      KIND = 4
      CALL FORMIC (KIND,ISN,JSN,IAUX,IC,LENGIC,B)

      DO 1 I = 1,LENGIC
        IF ( GETICM( IC(I),J,ICM,NICM ) ) THEN
          IC(I) = N2 + J
        ELSE
          WRITE (6,666) ISN,JSN,ICM,IC
  666     FORMAT ('0ERROR IN GETGIC--',17I4)
          CALL ABORT2
        ENDIF
    1 CONTINUE

      RETURN
      END
      SUBROUTINE NEWICM (NICMS)

*** INITIALIZE THE TABLE TO ZERO LENGTH

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)

      NICMS = 0

      RETURN
      END
      LOGICAL FUNCTION GETICM (ICM,IICM,ICMS,NICMS)

*** TABLE LOOKUP

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION ICMS(*)

      IICM = 1
  100 IF (IICM.LE.NICMS) THEN
        IF ( ICM.EQ.ICMS(IICM) ) THEN
          GETICM = .TRUE.
          RETURN
        ELSE
          IICM = IICM + 1
        ENDIF
        GOTO 100
      ENDIF

*** FELL THRU TABLE

      GETICM = .FALSE.

      RETURN
      END
      SUBROUTINE PUTICM (IC,L,ICMS,NICMS)

*** ADD ENTRIES OF IC INTO TABLE

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL GETICM
      DIMENSION ICMS(*)
      DIMENSION IC(31)

      MAXICM = (NVECS+1)*15+1

      DO 1 I = 1,L
        ICM = IC(I)
        IF ( GETICM(ICM,IICM,ICMS,NICMS) ) THEN
          CONTINUE
        ELSE
          IF (IICM.GT.MAXICM) THEN
            WRITE (6,666)
  666       FORMAT ('0ERROR - GPS INDEX OVERFLOW IN PUTICM')
            CALL ABORT2
          ELSE
            NICMS = IICM
            ICMS(NICMS) = ICM
          ENDIF
        ENDIF
    1 CONTINUE

      RETURN
      END
      SUBROUTINE RGPS (IUO,IOBS,NVEC,NR,NC,LENG,G,B,IAUX,IVF,SIGUWT,
     &                 A,NX,ITIME)

*** LIST ADJUSTED OBSERVATIONS AND RESIDUALS

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( MXSSN = 9999, NVECS = 16 )
      LOGICAL LBB,LGF,LCS,LVD,LVA,LVZ,LVS,
     &        LVR,LVG,LVC,LIS,LPS,LPG,LDR,LOS,LAP
      LOGICAL FATAL,NOBIGV
      LOGICAL LMSL,LSS,LUP
      CHARACTER*30 NAMES,NAME1,NAME2
      DIMENSION G(NR,NC),B(*),A(*),NX(*)
      DIMENSION ICM((NVECS+1)*15+1),GMDES(3*NVECS)
      DIMENSION KINDS(3*NVECS),ISNS(3*NVECS),JSNS(3*NVECS),LOBS(3*NVECS)
      COMMON /NAMTAB/ NAMES(MXSSN)
      COMMON /OPT/ AX,E2,DMSL,DGH,VM,VP,CTOL,ITMAX,ITMIN,IMODE,
     &             LMSL,LSS,LUP
      COMMON /OPRINT/ CRIT,LBB,LGF,LCS,LVD,LVA,LVZ,LVS,
     &                LVR,LVG,LVC,LIS,LPS,LPG,LDR,LOS,LAP
      SAVE ISNX,JSNX,ISN,JSN
      SAVE X1,Y1,Z1,X2,Y2,Z2,RNSUM
      SAVE NAME1,NAME2
      SAVE ITIMEX,JTIMEX
      DATA ISNX /0/, JSNX /0/, ITIMEX/0/, JTIMEX/0/

*** DO NOT ABORT/PRINT BIG RESIDUALS

      NOBIGV = .TRUE.

      READ (IUO,END = 666) ICM,NICM,KINDS,ISNS,JSNS,LOBS

*** LOAD WORK SPACE (G)

      DO 1 I = 1,NR
        READ (IUO,END=666) (G(I,J),J=1,NC)
    1 CONTINUE

*** COMPUTE AND DECORRELATE OBS. EQ. AND MISCLOSURE
*** IF MODE 3, COMPUTE STD. DEV. OF RESIDUALS AND REDUNDANCY NUMBERS
*** IF MODE 0, COMPUTE MARGINALLY DETECT. ERROR AND REDUNDANCY NUMBERS
*** IF MODE 3, COMPUTE MDE, PATCH WORK

      IF (IMODE.EQ.3 .OR. IMODE.EQ.0) THEN
        CALL GPSSDV (G,NR,NC,NVEC,LENG,B,ICM,NICM,KINDS,ISNS,JSNS,
     &        LOBS,IAUX,IVF,FATAL,NOBIGV,A,NX,N2,N3,N4,GMDES,ITIME)
      ELSE
        CALL GOBSEQ (G,NR,NC,NVEC,LENG,B,ICM,NICM,KINDS,ISNS,JSNS,
     &               LOBS,IAUX,IVF,FATAL,NOBIGV,N4,ITIME,N5)
      ENDIF

*** UNLOAD WORK SPACE

      IF (LVG) THEN
        CALL LINE2 (NR + 1)
        WRITE (6,2)
    2   FORMAT (' ')
      ENDIF

*** UNLOAD WORK SPACE

      RNSUM = 0.D0
      DO 500 I = 1,NR
        OBSB = G(I,NC)
        IF (IMODE.EQ.3) THEN
          VSD = G(I,N2)
          SD3 = G(I,N4)
          IF (SD3.EQ. 0.D0 ) SD3 = 1.D0 + 20
          RN = G(I,N3)
        ELSE
          VSD = G(I,N4)
          VSD3 = G(I,N4)
          RN = 0.D0
        ENDIF

        IF (IMODE.EQ.0 .OR. IMODE.EQ.3) THEN
          RN = G(I,N3)
          IF (RN.EQ. 0.D0 ) THEN
            GMDE = 0.D0
          ELSE
            GMDE = 3.D0 * DSQRT( GMDES(I) )
          ENDIF
        ENDIF

        KIND = KINDS(I)
        ISN = ISNS(I)
        JSN = JSNS(I)
        IOBS = LOBS(I)
        IF (ISN.NE.ISNX .OR. ITIME.NE.ITIMEX ) THEN
          ISNX = ISN
          ITIMEX = ITIME
          CALL GETXYZ (X1,Y1,Z1,B,ISN,ITIME)
          NAME1 = NAMES(ISN)
        ENDIF
        IF (JSN.NE.JSNX .OR. ITIME.NE.JTIMEX ) THEN
          JSNX = JSN
          JTIMEX = ITIME
          CALL GETXYZ (X2,Y2,Z2,B,JSN,ITIME)
          NAME2 = NAMES(JSN)
        ENDIF
        RNSUM = RNSUM + RN

*** GPS DEL X, KIND = 4

        IF (KIND.EQ.4) THEN
          OBS0 = X2 - X1
          IF (IAUX.GT.0) THEN
            CALL GETAUX (VAL,IAUX,B)
            OBS0 = OBS0 - OBS0 * VAL
          ENDIF
	  RX=OBS0
          IF (IMODE.EQ.0) OBSB = OBS0
          VX = OBS0 - OBSB
	  RBX=OBSB
          IF (IMODE.EQ.3) VSD3 = DIVID(VX,SD3)
          IF (SD3.EQ. 1.0D+20 ) SD3 = 0.D0
          CALL RSTAT (VSD,RN,VSD3,IOBS,KIND)
          IF (LSS) VSD = VSD / SIGUWT
	  IF (IMODE .EQ. 3) THEN
	    VSDX=SD3
	  ELSE
	    VSDX=VSD
	  ENDIF
          IF (LVG) THEN
            IF ( DABS(VSD) .GE.CRIT) THEN
              IF (IMODE.EQ.0) THEN
                WRITE (6,161) IOBS,OBS0,GMDE,NAME1,NAME2
  161           FORMAT (' ',I5,' DX',F15.4,F12.4,18X,
     &                  2(1X,A30) )
              ELSEIF (IMODE.EQ.3) THEN
                WRITE (6,162) IOBS,OBS0,OBSB,VX,SD3,VSD3,NAME1
  162           FORMAT (' ',I5,' DX',F15.4,F17.4,F10.4,F7.4,' METER',
     &                  F8.2,8X,A30)
              ELSE
                WRITE (6,16) IOBS,OBS0,OBSB,VX,VSD,NAME1,NAME2
   16           FORMAT (' ',I5,' DX',F15.4,F18.4,F10.4,
     &                  8X,' METER',F8.2,2(1X,A30) )
              ENDIF
            ENDIF
          ENDIF

*** DEL Y GPS, KIND = 5

        ELSEIF (KIND.EQ.5) THEN
          OBS0 = Y2 - Y1
          IF (IAUX.GT.0) THEN
            CALL GETAUX (VAL,IAUX,B)
            OBS0 = OBS0 - OBS0 * VAL
          ENDIF
	  RY=OBS0
          IF (IMODE.EQ.0) OBSB = OBS0
          VY = OBS0 - OBSB
	  RBY=OBSB
          IF (IMODE.EQ.3) VSD3 = DIVID(VY,SD3)
          IF (SD3.EQ. 1.0D+20 ) SD3 = 0.D0
          CALL RSTAT (VSD,RN,VSD3,IOBS,KIND)
          IF (LSS) VSD = VSD / SIGUWT
          IF (IMODE .EQ. 3) THEN 
            VSDY=SD3
          ELSE 
            VSDY=VSD 
          ENDIF
          IF (LVG) THEN
            IF ( DABS(VSD) .GE.CRIT) THEN
              IF (IMODE.EQ.0) THEN
                WRITE (6,171) IOBS,OBS0,GMDE
  171           FORMAT (' ',I5,' DY',F15.4,F12.4)
              ELSEIF (IMODE.EQ.3) THEN
                WRITE (6,172) IOBS,OBS0,OBSB,VY,SD3,VSD3,NAME2
  172           FORMAT (' ',I5,' DY',F15.4,F17.4,F10.4,F7.4,' METER',
     &                  F8.2,10X,A30)
              ELSE
                WRITE (6,17) IOBS,OBS0,OBSB,VY,VSD
   17           FORMAT (' ',I5,' DY',F15.4,F18.4,F10.4,8X,' METER'
     &                  ,F8.2)
              ENDIF
            ENDIF
          ENDIF

*** DEL Z GPS, KIND = 6

        ELSEIF (KIND.EQ.6) THEN
          OBS0 = Z2 - Z1
          IF (IAUX.GT.0) THEN
            CALL GETAUX (VAL,IAUX,B)
            OBS0 = OBS0 - OBS0 * VAL
          ENDIF
	  RZ=OBS0
          IF (IMODE.EQ.0) OBSB = OBS0
          VZ = OBS0 - OBSB
	  RBZ=OBSB
          IF (IMODE.EQ.3) VSD3 = DIVID(VZ,SD3)
          IF (SD3.EQ. 1.0D+20 ) SD3 = 0.D0
          CALL RSTAT (VSD,RN,VSD3,IOBS,KIND)
          IF (I.EQ.NR) THEN
            IF ( IMODE.EQ.3 .OR. IMODE.EQ.0 ) THEN
              CALL RSTAT2 (RNSUM)
            ENDIF
          ENDIF
          IF (LSS) VSD = VSD / SIGUWT
          IF (IMODE .EQ. 3) THEN 
            VSDZ=SD3
          ELSE 
            VSDZ=VSD 
          ENDIF
          CALL MNTMD(ITIME, IYR,IMON,IDAY,IHR,IMN)
	  IF (LVG) THEN
            IF ( DABS(VSD) .GE.CRIT) THEN
              IF (I.EQ.NR) THEN
                IF (IMODE.EQ.0) THEN
                  WRITE (6,181) IOBS,OBS0,GMDE,RNSUM,IHR,IMN,
     &                          IMON,IDAY,IYR
  181             FORMAT (' ',I5,' DZ',F15.4,F12.4,F18.2,2X,I2,':',
     &                     I2,1X,I2,'/',I2,'/',I4)
                ELSEIF (IMODE.EQ.3) THEN
                  WRITE (6,182) IOBS,OBS0,OBSB,VZ,SD3,VSD3,RNSUM,
     &               IHR,IMN,IMON,IDAY,IYR
  182             FORMAT (' ',I5,' DZ',F15.4,F17.4,F10.4,F7.4,' METER',
     &                    F8.2,F6.2,2X,I2,':',I2,1X,I2,'/',I2,'/',I4)
                ELSE
                  WRITE (6,18) IOBS,OBS0,OBSB,VZ,VSD
   18             FORMAT (' ',I5,' DZ',F15.4,F18.4,F10.4,8X,' METER',
     &                    F8.2)
                ENDIF
              ELSE
                IF (IMODE.EQ.0) THEN
                  WRITE (6,183) IOBS,OBS0,GMDE
  183             FORMAT (' ',I5,' DZ',F15.4,F12.4)
                ELSEIF (IMODE.EQ.3) THEN
                  WRITE (6,184) IOBS,OBS0,OBSB,VZ,SD3,VSD3
  184            FORMAT (' ',I5,' DZ',F15.4,F17.4,F10.4,F7.4,' METER',
     &                   F8.2)
                ELSE
                  WRITE (6,185) IOBS,OBS0,OBSB,VZ,VSD
  185             FORMAT (' ',I5,' DZ',F15.4,F18.4,F10.4,8X,' METER',
     &                    F8.2)
                ENDIF
              ENDIF
            ENDIF
	  ENDIF

*** FIND AND PRINT N, E, U RESIDUALS
 
        CALL TOLGH (RX,RY,RZ,R0N,R0E,R0U,ISN,JSN,B)
        CALL TOLGH (RBX,RBY,RBZ,RBN,RBE,RBU,ISN,JSN,B)
        CALL TOLGH (VX,VY,VZ,VN,VE,VU,ISN,JSN,B)   
        CALL COVLGH (VSDN,VSDE,VSDU,ISN,JSN,B,I,G,NR,NC,LENG)
	RESN=VN/VSDN
	RESE=VE/VSDE
	RESU=VU/VSDU
        IF (IMODE .EQ. 0) THEN
	  WRITE (*,*) ' '
          WRITE (*,191) ' DN',R0N
          WRITE (*,191) ' DE',R0E
          WRITE (*,191) ' DU',R0U
191       FORMAT (' ',6X,A3,F15.4)
	  WRITE (*,*) '           ------------------ '
        ELSEIF (IMODE .EQ. 3) THEN
          WRITE (*,*) ' '
          WRITE (*,192) ' DN',R0N,RBN,VN,VSDN,RESN
          WRITE (*,192) ' DE',R0E,RBE,VE,VSDE,RESE
          WRITE (*,192) ' DU',R0U,RBU,VU,VSDU,RESU
192       FORMAT (' ',5X,A3,F15.4,F17.4,F10.4,F7.4,' METER',F8.2)
	  WRITE (*,*) '           ------------------ '
        ELSE
          WRITE (*,*) ' '
          WRITE (*,193) ' DN',R0N,RBN,VN,RESN 
          WRITE (*,193) ' DE',R0E,RBE,VE,RESE
          WRITE (*,193) ' DU',R0U,RBU,VU,RESU
193       FORMAT (' ',6X,A3,F15.4,F18.4,F10.4,8x,' METER',  F8.2)
	  WRITE (*,*) '           ------------------ '
        ENDIF

*** ACCUMULATE N, E, U STATS, BYPASS DOWNWEIGHTED OUTLIERS AND NO CHECK
 
***        IF (DABS(VX+VY+VZ) .LE. 1.D-6) THEN
***          CONTINUE
***        ELSEIF (IMODE .EQ. 3 .AND. VSDZ .GT. 2.D0) THEN
***          CONTINUE
***        ELSEIF ((IMODE .EQ. 1 .OR. IMODE .EQ. 2) .AND. 
***     &       DABS(VSDZ) .LE. 0.02) THEN
***          CONTINUE
***        ELSE
***          CALL RSTAT4 (VN,VE,VU)
***        ENDIF
       ENDIF

  500 CONTINUE

      RETURN

*** PREMATURE END OF FILE

  666 WRITE (6,667) NVEC
  667 FORMAT ('0ERROR - PREMATURE FILE END IN RGPS -- NVEC=',I5)
      CALL ABORT2
      RETURN
      END
      SUBROUTINE GPSSDV (G,NR,NC,NVEC,LENG,B,ICM,NICM,KINDS,ISNS,JSNS,
     &        LOBS,IAUX,IVF,FATAL,NOBIGV,A,NX,N2,N3,N4,GMDES,ITIME)

*** COMPUTE AND DECORRELATE OBS. EQ. AND MISCLOSURE

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      PARAMETER ( NVECS = 16 )
      LOGICAL FATAL,NOBIGV,PROP,PROPCV
      LOGICAL LMSL,LSS,LUP
      DIMENSION KINDS(3*NVECS),ISNS(3*NVECS),JSNS(3*NVECS),LOBS(3*NVECS)
      DIMENSION ICM((NVECS+1)*15+1),CM((NVECS+1)*15+1)
      DIMENSION CM2((NVECS+1)*15+1),GMDES(3*NVECS)
      DIMENSION B(*),G(NR,NC),A(*),NX(*)
      DIMENSION IC(31),C(31)
      COMMON /NUMVFS/ NVFTOT,NVFREE
      COMMON /VFTB2/ VFS(30),VTV(30),VFRN(30),VFSS(30)
      COMMON /OPT/ AX,E2,DMSL,DGH,VM,VP,CTOL,ITMAX,ITMIN,IMODE,
     &             LMSL,LSS,LUP
      COMMON /LASTSN/ ISNX,JSNX,KSNX,ITIMEX,JTIMEX,KTIMEX
      COMMON /LASTXT/ ISNZ,JSNZ,KSNZ,ITIMEZ,JTIMEZ,KTIMEZ
      COMMON /STRUCT/ NSTA,NAUX,NUNK,IDIM,NSTAS,NOBS,NCON,NZ,NCD

*** LENGTH OF COEFF. ARRAYS

      LMAX = 6
      IF (NCD.GT.0) LMAX = LMAX + 24
      IF (IAUX.GT.0) LMAX = LMAX + 1

*** SET LAST STATION ENCOUNTERED TO NULL

      ISNX = 0
      JSNX = 0
      KSNX = 0
      ITIMEX = 0
      JTIMEX = 0
      KTIMEX = 0
      ISNZ = 0
      JSNZ = 0
      KSNZ = 0
      ITIMEZ = 0
      JTIMEZ = 0
      KTIMEZ = 0

*** DETERMINE WORK ARRAY COLUMN POINTERS

      CALL GLOCAT (N1,N2,N3,N4,N5,NVEC,LENG)

      NB = NUNK + 1
      NL = NR * (NR + 1) / 2
      NB1 = NX(NB)
      NB2 = NB1 + NL
      NB3 = NB2 + NL

*** (STEP 1.)  COMPUTE COEFFICIENTS

      DO 10 I = 1,NR
        LENGC = LMAX
        KIND = KINDS(I)
        ISN = ISNS(I)
        JSN = JSNS(I)
        CALL FORMC (KIND,C,B,ISN,JSN,IAUX,ITIME)
        CALL GETGIC (ISN,JSN,IAUX,ICM,IC,N2,NICM,B)
        CALL COMPIC(IC,C,LENGC)

        DO 2 J = 1,LENG
          K = J + N2
          G(I,K) = 0.D0
    2   CONTINUE

        DO 20 J = 1,LENGC
          G(I,IC(J) ) = C(J)
   20   CONTINUE
   10 CONTINUE

*** (STEP 2.)  PROPAGATE FOR STD. DEV. OF RESIDUALS
***            SCALE BY VARIANCE FACTOR AS NEEDED
***           (NOTE: COV. OF PARMS IS SCALED BY FACTORS IN SUB. NORMAL)

      DO 60 IROW = 1,NR
        DO 50 I = 1,NICM
          CM(I) = G(IROW,N2+I)
   50   CONTINUE

        IF ( .NOT.PROP(CM,ICM,NICM,VARL0,A,NX,IFLAG) ) THEN
          IF (IFLAG.EQ.1) THEN
            WRITE (6,666)
  666       FORMAT ('0STATE ERROR IN FIRST PROP OF GPSSDV')
            CALL ABORT2
          ELSE
            WRITE (6,667)
  667       FORMAT ('0PROFILE ERROR IN FIRST PROP OF GPSSDV')
            CALL ABORT2
          ENDIF
        ENDIF

        IDEX = INX(IROW,NR)
        A(IDEX+NB1) = VARL0
        IF (IROW.EQ.1) THEN
          VARLB = G(1,1) * G(1,1)
        ELSE
          VARLB = G(IROW,1)
        ENDIF

        IF (NVFTOT.GT.0) THEN
          IF (IVF.LT.0 .OR. IVF.GT.NVFTOT) THEN
            WRITE (6,663) IVF,NVFTOT
  663       FORMAT ('0ERROR - ILLEGAL IVF=',I10,' FOR N=',I10,
     &              ' IN GPSSDV')
            CALL ABORT2
          ELSEIF (IVF.NE.0) THEN
            VFACTR = VFS(IVF)
          ELSE
            VFACTR = 1.D0
          ENDIF
        ELSE
          VFACTR = 1.D0
        ENDIF
        VARLB = VARLB * VFACTR

        DIFF = VARLB - VARL0
        IF (DIFF.LT. 1.D-10 ) DIFF = 0.D0
        SDV = DSQRT(DIFF)
        G(IROW,N4) = SDV
   60 CONTINUE

*** (STEP 3.1)  LOAD ADJ. OBS. COV. MATRIX INTO 1ST SCRATCH SPACE

      DO 30 I = 1,NR
        DO 31 K = 1,NICM
          CM(K) = G(I,N2+K)
   31   CONTINUE
        IDEX = INX(I,NR) - I
        DO 33 J = I,NR
          IF (I.NE.J) THEN
            DO 32 K = 1,NICM
              CM2(K) = G(J,N2+K)
   32       CONTINUE
            IF ( .NOT.PROPCV(CM,ICM,NICM,CM2,ICM,NICM,
     &                                  COV,A,NX,IFLAG) ) THEN
              IF (IFLAG.EQ.1) THEN
                WRITE (6,668)
  668           FORMAT ('0STATE ERROR IN PROPCV OF GPSSDV')
                CALL ABORT2
              ELSE
                WRITE (6,669)
  669           FORMAT ('0PROFILE ERROR IN PROPCV OF GPSSDV')
                CALL ABORT2
              ENDIF
            ENDIF
            A(NB1+IDEX+J) = COV
          ENDIF
   33   CONTINUE
   30 CONTINUE

*** (STEP 3.2)  COPY OBS. COV. CHOLESKY FACTOR INTO 2ND SCRATCH SPACE

      IPOINT = NB2
      DO 40 I = 1,NR
        DO 41 J = I,NR
          IPOINT = IPOINT + 1
          A(IPOINT) = G(I,J)
   41   CONTINUE
   40 CONTINUE

*** (STEP 3.3)  INVERT IMAGE OF OBS. COV. CHOLESKY (A(NB3) IS SCRATCH)

      CALL AVERT(A(NB2+1),A(NB3+1),NR)

*** (STEP 4)  COMPUTE MARGINALLY DETECTABLE ERRORS (1 SIGMA)
***           TABLE GMDES(*) IS PART OF PATCH WORK

      IF (IMODE.EQ.0 .OR. IMODE.EQ.3) THEN
        DO 70 IROW = 1,NR
          DO 75 JCOL = 1,NICM
            CM(JCOL) = 0.D0
            DO 76 K = 1,NR
              IF (IROW.LE.K) THEN
                INDEX = INX(IROW,NR) - IROW + K
              ELSE
                INDEX = INX(K,NR) - K + IROW
              ENDIF
              CM(JCOL) = CM(JCOL) + G(K,N2+JCOL) * A(NB2+INDEX)
   76       CONTINUE
   75     CONTINUE
          IF ( .NOT.PROP(CM,ICM,NICM,VAR,A,NX,IFLAG) ) THEN
            IF (IFLAG.EQ.1) THEN
              WRITE (6,664)
  664         FORMAT ('0STATE ERROR IN SECOND PROP OF GPSSDV')
              CALL ABORT2
            ELSE
              WRITE (6,665)
  665         FORMAT ('0PROFILE ERROR IN SECOND PROP OF GPSSDV')
              CALL ABORT2
            ENDIF
          ENDIF
          IF (NVFTOT.GT.0) THEN
            IF (IVF.LT.0 .OR. IVF.GT.NVFTOT) THEN
              WRITE (6,663) IVF,NVFTOT
              CALL ABORT2
            ELSEIF (IVF.NE.0) THEN
              VFACTR = VFS(IVF)
            ELSE
              VFACTR = 1.D0
            ENDIF
          ELSE
            VFACTR = 1.D0
          ENDIF
          DEINV = A( NB2+INX(IROW,NR) ) * VFACTR - VAR
          IF (DEINV.LE. 1.D-6 ) THEN
            GMDES(IROW) = 0.D0
          ELSE
            GMDES(IROW) = 1.D0 / DEINV
          ENDIF
   70   CONTINUE
      ENDIF

*** (STEP 5.)  COMPUTE OBSERVATIONS AND MISCLOSURES

      IF (IMODE.NE.0) THEN
        DO 90 I = 1,NR
          KIND = KINDS(I)
          ISN = ISNS(I)
          JSN = JSNS(I)
          CALL COMPOB (KIND,OBS0,B,RDUMMY,ISN,JSN,IAUX,ITIME)
          CMO = OBS0 - G(I,NC)
          G(I,N2) = CMO
   90   CONTINUE
      ENDIF

*** (STEP 6.)  DECORRELATE MISCLOSURE

      IF (IMODE.NE.0) THEN
        CALL COMRHS (G,NR)

        IF (NVFTOT.GT.0 .AND. IVF.GT.0) THEN
          SFACTR = DSQRT(VFACTR)
          DO 15 I = 1,NR
            G(I,N2) = G(I,N2) / SFACTR
   15     CONTINUE
        ENDIF

*** TEST FOR LARGE MISCLOSURES

        IF (.NOT.NOBIGV) THEN
          SD = 1.D0
          DO 16 I = 1,NR
            CMO = G(I,N2)
            CALL BIGL (KINDS(I),ISNS(I),JSNS(I),LOBS(I),IVF,CMO,SD,
     &                 FATAL)
   16     CONTINUE
        ENDIF
      ENDIF

*** (STEP 7)  COMPUTE REDUNDANCY NUMBERS USING IMAGES
***           CAUTION--THIS STEP DAMAGES THE COEFFICIENTS IN G

      DO 80 I = 1,NR
        W = 0.D0
        DO 85 J = 1,NR
          IF (I.LE.J) THEN
            K = INX(I,NR) - I + J
          ELSE
            K = INX(J,NR) + I - J
          ENDIF
          SIGA = A(NB1+K)
          SIGL = A(NB2+K)
          W = W + SIGA * SIGL
   85   CONTINUE
        W = W / VFACTR
        RN = 1.D0 - W
        G(I,N3) = RN
   80 CONTINUE

      RETURN
      END
      SUBROUTINE AVERT (A,SCR,N)

*** INVERT TRIANGULAR CHOLESKY FACTOR IN PLACE

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DIMENSION A(*),SCR(*)

*** PUT CHOLESKY RECIPROCALS ON DIAGONAL

      DO 1 I = 1,N
        IDEX = INX(I,N)
        A(IDEX) = 1.D0 / A(IDEX)
    1 CONTINUE

*** COMPUTE GAUSSIAN FACTOR

      DO 2 I = 1,N
        DIAG = A( INX(I,N) )
        DO 6 J = I,N
          INDEX = INX(I,N) + J - I
          A(INDEX) = DIAG * A(INDEX)
    6   CONTINUE
    2 CONTINUE

*** PROCEED BY ROWS -- BOTTOM TO TOP

      DO 5 I = N-1, 1, -1
        I1 = I + 1
        IDEX = INX(I,N)
        DO 3 J = I1,N
          SCR(J) = -A(INX(I,N)+J-I)
    3   CONTINUE
        DO 4 J = I1,N
          INDEX = INX(I,N) + J - I
          A(INDEX) = 0.D0
          DO 8 K = I1,N
            IF (K.LE.J) THEN
              JNDEX = INX(K,N) + J - K
            ELSE
              JNDEX = INX(J,N) + K - J
            ENDIF
            A(INDEX) = A(INDEX) + SCR(K) * A(JNDEX)
    8     CONTINUE
    4   CONTINUE
        DO 7 J = I1,N
          A(IDEX) = A(IDEX) + SCR(J) * A(INX(I,N)+J-I)
    7   CONTINUE
    5 CONTINUE

      RETURN
      END
      INTEGER FUNCTION INX (I,N)

*** RETURN INDEX FOR DIAGONAL OF FULL TRIANGULAR MATRIX BY ROWS

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)

      NBAR = N * (N + 1) / 2
      J = N - I + 1
      JBAR = J * (J + 1) / 2
      INX = NBAR - JBAR + 1

      RETURN
      END
      SUBROUTINE RSTAT2 (RNSUM)

*** ACCUMULATE GPS NO-CHECK FOR MODE = 0 & MODE=3

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      COMMON /RESTAT/ VPV(32),RNUM(32),NUM(32),N0                   

      IF (RNSUM.LT. 3.0D-5 ) N0 = N0 + 3

      RETURN
      END
