      SUBROUTINE ADDAMB(CTIME,NPR,IPR,LOCPR,LOCSAT,LSATR,LSATU,
     2   DATAR,DATAU,UVECS,ZVEC,RNGMAX,LPLT,LBADP,
     3   NDATP,NAMB,IAMB,IAMBN,CAMBP)

C     ADDAMB uses a list of satellite pairs, along with information
C     about the locations of the satellites in the input structures,
C     and the locations of the pairs in the pairs database, to
C     calculate double-differenced ranges and L1 carrier phases.
C     It then calls AMBPROB using these double-differenced observables
C     as input.  Finally ADDAMB uses the results from AMBPROB to
C     increment the pairs database.

C     Arguments:
C        CTIME          [C*18]: Current epoch
C        NPR             [I*4]: Number of pairs to be processed.
C        IPR(2,*)        [I*4]: Array of the PRN numbers associated
C                               with the pairs.  The first entry is the
C                               lower PRN.
C        LOCPR(*)        [I*4]: Location of the i-th pair in the pairs
C                               database.
C        LOCSAT(32)      [I*4]: Location of PRN i in the arrays UVECS and
C                               ZVEC.
C        LSATR(*)        [I*4]: Location in DATAR of the satellite at
C                               location i in ISATC, UVECS, ZVEC, etc.
C        LSATU(*)        [I*4]: Location in DATAU of the satellite at
C                               location i in ISATC, UVECS, ZVEC, etc.
C        DATAR(*)        [C*?]: Structure containing the data for this
C                               epoch from the reference station.
C        DATAU(*)        [C*?]: Structure containint the data for this
C                               epoch from the reference station.
C        UVECS(3,*)      [R*8]: Array of unit vectors in the directions
C                               from the reference station to the various
C                               satellites at the current epoch.
C        ZVEC(3)         [R*8]: Unit vectors toward the zenith at the
C                               reference station
C        RNGMAX          [R*8]: Maximum possible double-differenced
C                               range [m].
C        LPLT            [L*1]: TRUE iff the user wants AMBPROB to
C                               generate an ambiguity-resolution plot for
C                               each pair.
C        LBADP(*)        [L*1]: Array element i set TRUE iff pair i is
C                               temporarily bad for this epoch.
C        NDATP(*)        [I*4]: Number of epochs processed for the pair
C                               in slot i of the pairs database.
C        NAMB(*)         [I*4]: Number of different ambiguities found so
C                               far for the pair in slot i of the pairs
C                               database.
C        IAMB(MAXAMB,*)  [I*4]: Values of the double-differenced
C                               integer-cycle corrections found so far
C                               for the pair in slot i of the pairs
C                               database.
C        IAMBN(MAXAMB,*) [I*4]: The number of times the integer-cycle
C                               correction in the corresponding location
C                               in IAMB has be found to be most probable.
C        CAMBP(MAXAMB,*) [R*8]: Summed probability of the integer-cycle
C                               correction at the corresponding location
C                               in IAMB.


      IMPLICIT REAL*8 (A-H,O-Z)
      LOGICAL*1 LBADP,LPLT
      CHARACTER CTIME*18

      PARAMETER (MAXSAT = 10, MAXAMB = 50, MXXSAT = 32,
     2   C = 2.99792458D8, FREQ1 = 1.57542D9, FREQ2 = 1.22760D9)
      PARAMETER (WVLTH1 = C/FREQ1, WVLTH2 = C/FREQ2)

      DIMENSION CAMBP(MAXAMB,*),CPHASE(2),CPERR(2),CPROB(20),
     2   CRANGE(20),DATAR(MAXSAT+1),DATAU(MAXSAT+1),IAMB(MAXAMB,*),
     3   IAMBN(MAXAMB,*),IPR(2,*),LBADP(*),LOCPR(*),LOCSAT(MXXSAT),
     4   LSATR(*),LSATU(*),NAMB(*),NDATP(*),PRANGE(2),PRERR(2),
     5   UVECS(3,*),VARMPC(2),VARMPR(2),ZMPC(2),ZMPR(2),ZVEC(3)

      STRUCTURE /CHDATA/         ! Receiver data for each channel (L1, L2)
         REAL*4 SNR              ! Voltage SNR
         REAL*8 DELAY            ! Pseudo-range [sec]
         REAL*4 SIGDLA           ! Uncertainty of pseudo-range [sec]
         REAL*8 PHASE            ! Phase [cycles]
         REAL*4 SIGPHS           ! Uncertainty of phase [cycles]
      END STRUCTURE

      STRUCTURE /ROGDATA/        ! Full data record for 1 satellite, 1 epoch
         INTEGER*4 ISEC50,ICHAN,ISAT,INT,IRSTAT ! UT at center of interval,
         RECORD /CHDATA/ P1, P2           ! channel #, PRN #, integration
      END STRUCTURE                       ! time [sec], status, P1 and P2 data

      RECORD /ROGDATA/ DATAR,DATAU        ! Reference and User stations

      DATA CPERRLIM /6.0D-2/  ! Largest allowed double-diff phase error [cyc.]
      DATA ZMPR /2.5D-1,1.0D+1/        ! Zenith multipath for L1, L2 range [m]
      DATA ZMPC /1.0D-3,2.0D-3/        ! Zenith multipath for L1, L2 phase [m]


      DO  JPR=1,NPR                                   ! For each pair
         IF (.NOT.LBADP(JPR))  THEN
            LS1 = LOCSAT(IPR(1,JPR))               ! Locations of the two
            LS2 = LOCSAT(IPR(2,JPR))               ! satellites in ISATC
            LSR1 = LSATR(LS1)                      ! Locations of satellites
            LSR2 = LSATR(LS2)                      ! 1 and 2 of the pair in
            LSU1 = LSATU(LS1)                      ! DATAR and DATAU
            LSU2 = LSATU(LS2)
               ! COMPUTE MULTIPATH CONTRIBUTION TO DOUBLE-DIFF. VARIANCE
            SINEL1 = UVECS(1,LS1)*ZVEC(1) + UVECS(2,LS1)*ZVEC(2) +
     2               UVECS(3,LS1)*ZVEC(3)       ! Sin(elev.) for first sat.
            SINEL2 = UVECS(1,LS2)*ZVEC(1) + UVECS(2,LS2)*ZVEC(2) +
     2               UVECS(3,LS2)*ZVEC(3)       ! Sin(elev.) for second sat.
            VARMPR(1) = 2.0D0*((ZMPR(1)/SINEL1)**2 +  ! Multipath variance for
     2                         (ZMPR(1)/SINEL2)**2)/(C*C)   ! L1 range [s^2]
            VARMPR(2) = ((ZMPR(2)/ZMPR(1))**2)*VARMPR(1)    ! Also L2
            VARMPC(1) = 2.0D0*((ZMPC(1)/SINEL1)**2 +  ! Same for carrier phase
     2                         (ZMPC(1)/SINEL2)**2)/(WVLTH1*WVLTH1)  ! [cyc.]
            VARMPC(2) = ((ZMPC(2)/ZMPC(1))**2)*VARMPC(1)    ! And L2
            PRANGE(1) = C*((DATAU(LSU2).P1.DELAY -       ! Form double-diff.
     2                      DATAR(LSR2).P1.DELAY) -      ! of P1 pseudo-range
     3                     (DATAU(LSU1).P1.DELAY -       ! and convert from
     4                      DATAR(LSR1).P1.DELAY))       ! seconds to meters
            PRANGE(2) = C*((DATAU(LSU2).P2.DELAY -
     2                      DATAR(LSR2).P2.DELAY) -      ! Same for P2
     3                     (DATAU(LSU1).P2.DELAY -
     4                      DATAR(LSR1).P2.DELAY))
            PRERR(1) = C*SQRT(DATAU(LSU2).P1.SIGDLA**2 +       ! Compute
     2                        DATAR(LSR2).P1.SIGDLA**2 +       ! errors in
     3                        DATAU(LSU1).P1.SIGDLA**2 +       ! corresponding
     4                        DATAR(LSR1).P1.SIGDLA**2 + VARMPR(1))
            PRERR(2) = C*SQRT(DATAU(LSU2).P2.SIGDLA**2 +       ! elements of
     2                        DATAR(LSR2).P2.SIGDLA**2 +       ! PRANGE
     3                        DATAU(LSU1).P2.SIGDLA**2 +
     4                        DATAR(LSR1).P2.SIGDLA**2 + VARMPR(2))
            CPHASE(1) = (DATAU(LSU2).P1.PHASE -
     2                   DATAR(LSR2).P1.PHASE) -            ! Compute double-
     3                  (DATAU(LSU1).P1.PHASE -             ! differenced L1
     4                   DATAR(LSR1).P1.PHASE)              ! phase [cycles]
            CPHASE(2) = (DATAU(LSU2).P2.PHASE -
     2                   DATAR(LSR2).P2.PHASE) -
     3                  (DATAU(LSU1).P2.PHASE -             ! Same for L2
     4                   DATAR(LSR1).P2.PHASE)
            CPERR(1) = SQRT(DATAU(LSU2).P1.SIGPHS**2 +      ! Compute
     2                      DATAR(LSR2).P1.SIGPHS**2 +      ! errors in CPHASE
     3                      DATAU(LSU1).P1.SIGPHS**2 +
     4                      DATAR(LSR1).P1.SIGPHS**2 + VARMPC(1))
            IF (CPERR(1).GT.CPERRLIM)  THEN     ! Don't let huge value through
               WRITE(6,'(''*** PAIR'',I3,'' ('',I2.2,''-'',I2.2,
     2            ''): CPERR(1) ='',1PD11.4,'' CYCLES: REDUCE TO'',
     3            D11.4,'' ***'')' )
     4               JPR, IPR(1,JPR),IPR(2,JPR), CPERR(1),CPERRLIM
               CPERR(1) = CPERRLIM
            END IF
            CPERR(2) = SQRT(DATAU(LSU2).P2.SIGPHS**2 +
     2                      DATAR(LSR2).P2.SIGPHS**2 +
     3                      DATAU(LSU1).P2.SIGPHS**2 +
     4                      DATAR(LSR1).P2.SIGPHS**2 + VARMPC(2))
            IF (CPERR(2).GT.CPERRLIM)  THEN     ! Don't let huge value through
               WRITE(6,'(''*** PAIR'',I3,'' ('',I2.2,''-'',I2.2,
     2            ''): CPERR(2) ='',1PD11.4,'' CYCLES: REDUCE TO'',
     3            D11.4,'' ***'')' )
     4               JPR, IPR(1,JPR),IPR(2,JPR), CPERR(2),CPERRLIM
               CPERR(2) = CPERRLIM
            END IF
                           ! WRITE RESIDUALS TO AN OUTPUT FILE, IF DESIRED
            CALL AMBPROB(PRANGE,PRERR,CPHASE,CPERR,1.0D-1,RNGMAX,
     2                   LPLT, NOUT,CRANGE,CPROB,CRERR)     ! Find ambiguities
            LOCDB = LOCPR(JPR)      ! Location of this pair in pairs database
            NDATP(LOCDB) = NDATP(LOCDB) + 1     ! Incr. # of epochs this pair
            DO  I=1,NOUT               ! For each possible ambiguity returned
               JAMB = JIDNNT(CPHASE(1) - CRANGE(I)/WVLTH1)        ! Ambiguity
               J = 1                            ! Initialize search location
               DO  WHILE ((IAMB(J,LOCDB).NE.JAMB).AND.   ! Look through
     2                    (J.LE.NAMB(LOCDB)))         ! existing ambiguities
                  J = J + 1                           ! for the current value
               END DO
               IF (J.LE.NAMB(LOCDB))  THEN      ! Current ambiguity is in list
                  IF (I.EQ.1)  IAMBN(J,LOCDB) = IAMBN(J,LOCDB) + 1
                  CAMBP(J,LOCDB) = CAMBP(J,LOCDB) + CPROB(I)
               ELSE IF (NAMB(LOCDB).LT.MAXAMB)  THEN     ! New entry in list
                  NAMB(LOCDB) = J            ! (Which equals NAMB(LOCDB) + 1)
                  IAMB(J,LOCDB) = JAMB       ! Assign new ambiguity value
                  IF (I.EQ.1)  THEN          ! Initialize count of number of
                     IAMBN(J,LOCDB) = 1      ! times this ambiguity has ranked
                  ELSE                       ! #1, depending on whether it's
                     IAMBN(J,LOCDB) = 0      ! #1 this time
                  END IF
                  CAMBP(J,LOCDB) = CPROB(I)           ! Init. probability sum
               ELSE                 ! We need a new entry but there is no room
                  TYPE '('' * >> AMBIG. AT '',A,
     2               '', PR ('',I2.2,'','',I2.2,''): IAMB ='',
     3               I10,'', P ='',1PD9.2,'' *'')',         ! Warn & forget it
     4                  CTIME, IPR(1,JPR),IPR(2,JPR), JAMB, CPROB(I)
               END IF
            END DO
         END IF
      END DO

      RETURN
      END
