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