<<<<<<< sumlin.for
* Last processed by NICE on 18-Sep-1995 09:56:00                    
* Customized for :  EEEI        HPUX        UNIX         
=======
* Last processed by NICE on 24-Oct-1995 15:39:00                    
* Customized for :  EEEI        HPUX        UNIX         
>>>>>>> 1.16
      SUBROUTINE SUMLIN (ERROR,USER_FUNCTION)
C----------------------------------------------------------------------
C SAS	Support routine for command SUM
C	Add all scans of current index together with various
C	weights and options. Update the origin information and
C	reset the plotting constants
C Arguments :
C	ERROR	L	Logical error flag	Output
C (5-mar-1985)
C----------------------------------------------------------------------
      LOGICAL ERROR, FSEC, SIC_CTRLC, USER_FUNCTION
      EXTERNAL USER_FUNCTION
*
* Add all scans in current index together
      INCLUDE 'inc:memory.inc'
      INCLUDE 'common.inc'
      INCLUDE 'index.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'rdata.inc'
      INCLUDE 'structure.par'
      INCLUDE 'setup.inc'
*
<<<<<<< sumlin.for
      INTEGER KX,NDATA,NSUM,I,N, SIC_GETVM, SAVE_VTYPE
      INTEGER(KIND=4) AW1,AW2,ADDR,POINTER
      CHARACTER MESS*80, ANSWER*4, QUIT*4
      REAL RSUM
      LOGICAL SWAP, SWITCH, NOTOK
=======
      INTEGER KX,NDATA,NSUM,I,N, SIC_GETVM, SAVE_VTYPE, IR
      INTEGER(KIND=4) AW1,AW2,ADDR,POINTER
      CHARACTER NOMBRE*80, ANSWER*4, QUIT*4
      REAL RSUM, TIME_SAVED, TSYS_SAVED
      LOGICAL SWAP, SWITCH, NOTOK, END
>>>>>>> 1.16
* Data
      INCLUDE 'inc:constant.inc'
      DATA QUIT/'QUIT'/
*
* Verify if any input file
      IF (ILUN.EQ.0) THEN
         CALL MESSAGE(8,3,'SUM','No input file connected')
         ERROR = .TRUE.
         RETURN
      ENDIF
      IF (CXNEXT.LE.1) THEN
         CALL MESSAGE(8,3,'SUM','Index is empty')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Patch COMPOSITE
      IF (SALIG2.NE.'I') THEN
         ERROR = SIC_GETVM(2*MDATA,ADDR).NE.1
         IF (ERROR) RETURN
         AW1 = POINTER(ADDR,MEMORY)
         AW2 = AW1+MDATA
      ENDIF
*
      NSUM = 0
      KNEXT= 1
      KX   = CX_IND(1)
      CALL COPYRT(USER_FUNCTION)
*
* read the first observation into R
*
      CALL GET_FIRST(USER_FUNCTION,ERROR)
      TIME_SAVED = RTIME
      TSYS_SAVED = RTSYS
      NSUM = 1
<<<<<<< sumlin.for
      WRITE(MESS,31) KX,LNUM,ABS(LVER),LSCAN
      CALL MESSAGE(7,1,'SUM',MESS)
      CALL RZERO(USER_FUNCTION)    ! initialise le buffer R
      DO I=-20,0
         PRESEC(I) = FSEC(I)
      ENDDO
      CALL RGEN(ERROR)             ! infos generales
      ERROR=.FALSE.
      CALL RPOS(ERROR)             ! positions
      ERROR = .FALSE.
      CALL CONVERT_POS
      CALL RSPEC(ERROR)            ! spectro
      ERROR = .FALSE.
      CALL CONVERT_VEL
      IF (SBEAMT.NE.0) THEN
         CALL RCAL(ERROR)          ! Calibration section
         IF (ERROR) THEN
            MESS = 'Use SET CAL OFF to disable calibration '
     $      //'checks in averaging'
            CALL MESSAGE(7,2,'SUM',MESS)
            GOTO 20
         ENDIF
      ENDIF
      IF (PRESEC(BASE_SEC)) THEN
         CALL RBASE(ERROR)         ! ligne de base
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(PLOT_SEC)) THEN
         CALL RPLOT(ERROR)         ! limites trace
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(ORIG_SEC)) THEN
         CALL RORIG(ERROR)         ! origine (n0s de manips)
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(GAUS_SEC)) THEN
         CALL RGAUS(ERROR)         ! parametres fit gaussiennes
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(SHELL_SEC)) THEN
         CALL RSHEL(ERROR)         ! parametres fit gaussiennes
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(NH3_SEC)) THEN
         CALL RNH3(ERROR)          ! parametres fit gaussiennes
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(ABS_SEC)) THEN
         CALL RABS(ERROR)          ! parametres fit gaussiennes
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(SWITCH_SEC)) THEN
         CALL RFSW(ERROR)          ! parametres F-SWITCH
         ERROR = .FALSE.
      ELSE
         RSWMOD = MOD_POS
      ENDIF
=======
d      WRITE(NOMBRE,31) KX,LNUM,ABS(LVER),LSCAN
d      CALL MESSAGE(7,1,'SUM',NOMBRE)
*
>>>>>>> 1.16
      IF (RSWMOD.EQ.MOD_FREQ) THEN
         SWITCH = .TRUE.
      ELSE
         SWITCH = .FALSE.
      ENDIF
<<<<<<< sumlin.for
C+VMS+UNIX-SUNOS
      NOTOK = USER_FUNCTION('GET')
      IF (NOTOK) THEN
         CALL MESSAGE(4,2,'GET','Error reading user sections')
      ENDIF
C-VMS-UNIX-SUNOS
      NDATA=MDATA
      CALL RDATA(NDATA,RDATAY,ERROR)
      ERROR = .FALSE.
      RXNUM = KX
=======
>>>>>>> 1.16
*
* Patch COMPOSITE
      IF (SALIG2.NE.'I') THEN
         CALL ADDINIT (MEMORY(AW1),ERROR)
         SWAP = .FALSE.
      ENDIF
* Temporily change some default settings
      IF (SVELOC.EQ.VEL_AUT) THEN
         SAVE_VTYPE = SVELOC
         SVELOC = RVTYPE
      ELSE
         SAVE_VTYPE = VEL_LSR      ! Anything but VEL_AUT would be ok
      ENDIF
<<<<<<< sumlin.for
*-------
=======
>>>>>>> 1.16
*
* Loop here on next observation
      END = .FALSE.
      DO WHILE (.NOT.END)
         DO IR=1, R_NDUMP
*
            IF (NSUM.GT.1) THEN
               IF (RSWMOD.EQ.MOD_FREQ) THEN
                  IF (.NOT.SWITCH) CALL MESSAGE(6,2,'SUM',
     $            'Adding folded to unfolded spectra')
                  SWITCH = .TRUE.
               ELSEIF (SWITCH) THEN
                  CALL MESSAGE(6,2,'SUM',
     $            'Adding folded to unfolded spectra')
                  RNPHAS = 0
               ENDIF
               IF (SALIG2.EQ.'I') THEN
                  CALL ADDLIN (ERROR)
*
<<<<<<< sumlin.for
      CALL ROBS(KX,ERROR)          ! lit le header de l'observation
      IF (ERROR) GOTO 20
      NSUM = NSUM + 1
      WRITE(MESS,31) KX,LNUM,ABS(LVER),LSCAN
31    FORMAT('Entry ',I4,' Observation ',I5,';',I3,' Scan ',I5)
      CALL MESSAGE(7,1,'SUM',MESS)
      CALL RZERO(USER_FUNCTION)    ! initialise le buffer R
      DO I=-20,0
         PRESEC(I) = FSEC(I)
      ENDDO
      CALL RGEN(ERROR)             ! infos generales
      ERROR=.FALSE.
      CALL RPOS(ERROR)             ! positions
      ERROR = .FALSE.
      CALL CONVERT_POS
      CALL RSPEC(ERROR)            ! spectro
      ERROR = .FALSE.
      CALL CONVERT_VEL
      IF (SBEAMT.NE.0) THEN
         CALL RCAL(ERROR)          ! Calibration section
         IF (ERROR) GOTO 20
      ENDIF
      IF (PRESEC(BASE_SEC)) THEN
         CALL RBASE(ERROR)         ! ligne de base
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(PLOT_SEC)) THEN
         CALL RPLOT(ERROR)         ! limites trace
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(ORIG_SEC)) THEN
         CALL RORIG(ERROR)         ! origine (n0s de manips)
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(GAUS_SEC)) THEN
         CALL RGAUS(ERROR)         ! parametres fit gaussiennes
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(SHELL_SEC)) THEN
         CALL RSHEL(ERROR)         ! parametres fit gaussiennes
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(NH3_SEC)) THEN
         CALL RNH3(ERROR)          ! parametres fit gaussiennes
         ERROR = .FALSE.
      ENDIF
      IF (PRESEC(ABS_SEC)) THEN
         CALL RABS(ERROR)          ! parametres fit gaussiennes
         ERROR = .FALSE.
      ENDIF
* Switching parameters
      IF (PRESEC(SWITCH_SEC)) THEN
         CALL RFSW(ERROR)
         ERROR = .FALSE.
      ELSE
         RSWMOD = MOD_POS
      ENDIF
      IF (RSWMOD.EQ.MOD_FREQ) THEN
         IF (.NOT.SWITCH) CALL MESSAGE(6,2,'SUM',
     $   'Adding folded to unfolded spectra')
         SWITCH = .TRUE.
      ELSEIF (SWITCH) THEN
         CALL MESSAGE(6,2,'SUM','Adding folded to unfolded spectra')
         RNPHAS = 0
      ENDIF
*
C+VMS+UNIX-SUNOS
      NOTOK = USER_FUNCTION('GET')
      IF (NOTOK) THEN
         CALL MESSAGE(4,2,'GET','Error reading user sections')
      ENDIF
C-VMS-UNIX-SUNOS
      NDATA = MDATA
      CALL RDATA(NDATA,RDATAY,ERROR)
      ERROR = .FALSE.
      RXNUM = KX
*
      IF (SALIG2.EQ.'I') THEN
         CALL ADDLIN (ERROR)
*
=======
>>>>>>> 1.16
* Patch COMPOSITE
               ELSE
                  IF (SWAP) THEN
                     CALL ADDOVER (MEMORY(AW2),MEMORY(AW1),ERROR)
                  ELSE
                     CALL ADDOVER (MEMORY(AW1),MEMORY(AW2),ERROR)
                  ENDIF
                  SWAP = .NOT.SWAP
               ENDIF
               IF (ERROR) GOTO 20
            ENDIF
*
* get next record as a new spectrum:
            IF (IR.LT.R_NDUMP) THEN
               CALL COPYRT(USER_FUNCTION)
               CALL GET_REC(ERROR)
               IF (ERROR) GOTO 20
               R_RECORD = R_RECORD+1
               RTIME = TIME_SAVED
               RTSYS = TSYS_SAVED
               NSUM = NSUM+1
            ENDIF
         ENDDO
*
* Loop over index
         IF (SIC_CTRLC()) THEN
            DO WHILE (.TRUE.)
               CALL SIC_WPRN('W-SUM,  <^C> pressed, '//
     $         'type Q to abort, RETURN to continue'
     $         ,ANSWER,N)
               CALL SIC_BLANC(ANSWER,N)
               IF (N.EQ.0) GOTO 10
               N = MIN(N,4)
               IF (ANSWER(1:N).EQ.QUIT(1:N)) THEN
                  CALL MESSAGE(8,2,'SUM','Sum interrupted by <^C>')
                  ERROR = .TRUE.
                  RETURN
               ENDIF
            ENDDO
         ENDIF
10       CALL GET_NEXT(END,USER_FUNCTION,ERROR)
         IF (ERROR) GOTO 20
         TIME_SAVED = RTIME
         TSYS_SAVED = RTSYS
         NSUM = NSUM + 1
d         WRITE(NOMBRE,31) KX,LNUM,ABS(LVER),LSCAN
d 31      FORMAT('Entry ',I4,' Observation ',I5,';',I3,' Scan ',I5)
d         CALL MESSAGE(7,1,'SUM',NOMBRE)
      ENDDO
*
20    CONTINUE
*
* Patch COMPOSITE
      IF (SALIG2.NE.'I') THEN
         CALL FREE_VM(2*MDATA,ADDR)
*
      ELSEIF (NSUM.GT.1) THEN
         RSIGFI = 0.0
         IF (SWEIGH.EQ.'E') THEN
            RSUM = 1./NSUM
            DO I=1,RNCHAN
               IF (RDATAY(I).NE.RBAD) RDATAY(I) = RDATAY(I)*RSUM
            ENDDO
         ENDIF
      ENDIF
*
* The summed spectrum is generated by program:
      RXNUM = -1
*
* The summed spectrum has a single record:
      R_NDUMP = 0
      PRESEC(DESCR_SEC) = .FALSE.
      PRESEC (ORIG_SEC) = .TRUE.
      CALL NEWDAT
      IF (SAVE_VTYPE.EQ.VEL_AUT) THEN
         SVELOC = SAVE_VTYPE
      ENDIF
      END
*<FF>
      SUBROUTINE ADDLIN (ERROR)
C----------------------------------------------------------------------
C SAS	Internal routine
C	Adds buffer R and T together.
C	The position compatibility checked if SMATCH is true
C	alignment with respect to frequency if SALIGN = 'F'
C	                          velocity if SALIGN = 'V'
C				  channel number if SALIGN = 'C' (default)
C
C Arguments :
C	ERROR	L	Logical error flag	Output
C----------------------------------------------------------------------
      REAL*8 C
      PARAMETER (C=299792.458D0)
      INCLUDE 'index.inc'
      INCLUDE 'parameter.inc'
      INCLUDE 'rdata.inc'
      INCLUDE 'structure.par'
      INCLUDE 'tdata.inc'
      INCLUDE 'structure.pat'
      INCLUDE 'setup.inc'
      INCLUDE 'work.inc'
      LOGICAL ERROR
      REAL*4 S(MDATA)
      EQUIVALENCE (S,DATA_BUFFER)
      INTEGER*4 SNCHAN
*
* Variables locales
      REAL*4 DPOS2, STOL2, AT, AR, BT, BR, XT1, XT2, XR1, XR2,
     $XRMIN, XRMAX, XTMIN, XTMAX, XS1, XS2, TW, RW, SW,
     $AS, BS, RDXMAX, RDXMIN, TDXMAX, TDXMIN, XS, WW, SS,
     $XR, DX, W, XT, STIME, SVOFF, STSYS, DPOSR, DPOST
      INTEGER IS, IR, IT, I, ITMIN, ITMAX, IRMIN, IRMAX
      LOGICAL ALL
      CHARACTER*80 CHAIN, MESS
      INTEGER*4 SNSEQ
      INTEGER*4 SSTART(0:MSEQ),SEND(0:MSEQ)
*
* Kind compatibility
*
      IF (RKIND.NE.TKIND) THEN
         CALL MESSAGE(8,3,'ADD',
     $   'Cannot add Continuum with Line data')
         ERROR = .TRUE.
         RETURN
      ELSEIF (PRESEC(XCOORD_SEC).OR.TSECPR(XCOORD_SEC)) THEN
         MESS = 'Irregularly sampled data not yet supported'
         CALL MESSAGE(6,3,'ADD',MESS)
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Position compatibility
*
      IF (SMATCH) THEN
         STOL2 = STOLE ** 2
         IF (RTYPEC.NE.TTYPEC) THEN
            DPOSR = (RLAMOF*COS(RBET))**2+RBETOF**2
            DPOST = (TLAMOF*COS(TBET))**2+TBETOF**2
            IF (DPOSR.GT.STOL2 .OR. DPOST.GT.STOL2) THEN
               CALL MESSAGE(7,3,'ACCUMULATE',
     $         'Coordinate systems are not compatible')
               ERROR = .TRUE.
               RETURN
            ENDIF
         ENDIF
         DPOS2 = ((RLAM+RLAMOF-TLAM-TLAMOF)*COS(RBET))**2
     $   + (RBET+RBETOF-TBET-TBETOF)**2
         IF (DPOS2 .GT. STOL2) THEN
            CALL MESSAGE(7,3,'ACCUMULATE',
     $      'Positions are not compatible')
            ERROR = .TRUE.
            RETURN
         ENDIF
      ENDIF
*
* Calibration Compatibility
*
      IF (SBEAMT.NE.0) THEN
         IF (ABS(RBEEFF-TBEEFF).GT.SBEAMT) THEN
            WRITE(CHAIN,101) RBEEFF,TBEEFF
            MESS = 'Different calibrations'//CHAIN
            CALL MESSAGE(8,3,'ACCUMULATE',MESS)
            ERROR = .TRUE.
            RETURN
         ENDIF
         IF (SGAINT.NE.0) THEN
            IF (ABS(RGAINI-TGAINI).GT.SGAINT) THEN
               WRITE(CHAIN,101) RBEEFF,TBEEFF
               MESS = 'Different calibrations'//CHAIN
               CALL MESSAGE(8,3,'ACCUMULATE',MESS)
               ERROR = .TRUE.
               RETURN
            ENDIF
         ENDIF
      ENDIF
*
* Frequency compatibility
*
      IF (RVTYPE.NE.TVTYPE) THEN
         MESS = 'Different velocity definition for R and T'
         CALL MESSAGE(6,2,'ACCUMULATE',MESS)
      ENDIF
      IF (SALIGN.EQ.'C') THEN
         AT = 1.
         AR = 1.
         BT = TVOFF - TVRES*TRCHAN + (RRESTF-TRESTF)/RRESTF*C
         BR = RVOFF - RVRES*RRCHAN
         IF (RRCHAN.NE.TRCHAN) THEN
            WRITE (CHAIN,101) RRCHAN,TRCHAN
            MESS = 'Different reference channels'//CHAIN
            CALL MESSAGE(6,2,'ACCUMULATE',MESS)
         ENDIF
         IF (ABS(RFRES).NE.ABS(TFRES)) THEN
            WRITE(CHAIN,101) RFRES,TFRES
            MESS = 'Different resolutions'//CHAIN
            CALL MESSAGE(6,2,'ACCUMULATE',MESS)
         ENDIF
         IF (ABS(BT-BR).GT.ABS(RVRES)*1E-2) THEN     ! SG 22-Apr-1986
            CALL MESSAGE(6,2,'ACCUMULATE',
     $      'Spectra not aligned in sky frequency')
            WRITE(6,*) BR,BT
         ENDIF
         BT = 0.
         BR = 0.
      ELSEIF ( SALIGN.EQ.'F') THEN
         AT = TFRES
         AR = RFRES
         BR = - RRCHAN * RFRES
         BT = - TRCHAN * TFRES + (TRESTF - RRESTF)
         IF (RVOFF.NE.TVOFF) THEN
            WRITE(CHAIN,101) RVOFF,TVOFF
            MESS = 'Different velocities'//CHAIN
            CALL MESSAGE(6,2,'ACCUMULATE',MESS)
         ENDIF
      ELSEIF (SALIGN.EQ.'V') THEN
         AT = TVRES
         AR = RVRES
         BT = TVOFF - TVRES*TRCHAN
         BR = RVOFF - RVRES*RRCHAN
         IF (RRESTF.NE.TRESTF) THEN
            WRITE(CHAIN,101) RRESTF,TRESTF
            MESS = 'Different rest frequencies'//CHAIN
            CALL MESSAGE(6,2,'ACCUMULATE',MESS)
         ENDIF
      ENDIF
*
* Bad channel summation
      ALL = SBAD.EQ.'O'
*
      XT1 = AT + BT
      XT2 = AT*TNCHAN + BT
      XR1 = AR + BR
      XR2 = AR*RNCHAN + BR
      XTMIN = MIN(XT1,XT2)
      XTMAX = MAX(XT1,XT2)
      XRMIN = MIN(XR1,XR2)
      XRMAX = MAX(XR1,XR2)
*
* Cas Intersection
      IF (SALIG2.EQ.'I') THEN
         XS1 = MAX(XTMIN,XRMIN)
         XS2 = MIN(XTMAX,XRMAX)
         IF (XS1 .GE. XS2) THEN
            MESS = SALIGN//' ranges do not intersect'
            CALL MESSAGE(7,3,'ACCUMULATE',MESS)
            ERROR = .TRUE.
            RETURN
         ENDIF
      ELSE
* Cas reunion
         XS1 = MIN(XTMIN,XRMIN)
         XS2 = MAX(XTMAX,XRMAX)
      ENDIF
*
* Determine parameters of sum
*
      AS = MAX(ABS(AR),ABS(AT))
      SNCHAN = (XS2 - XS1)/AS + 1.5
      IF (SNCHAN.GT.MDATA) THEN
         CALL MESSAGE(8,4,'ACCUMULATE','Too many channels requested')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Preference is given to previously accumulated spectra if resolutions are
* identical, in order to avoid a slow erosion in the case of accumulation of
*  many spectra, with not precisely aligned channels.
*
      IF (ABS(AS-ABS(AT)).GT.1.E-6*AS) THEN
         XS1 = AR * NINT((XS1-BR)/AR) + BR
         SVOFF = RVOFF
      ELSE
         XS1 = AT * NINT((XS1-BT)/AT) + BT
         SVOFF = TVOFF
      ENDIF
      BS = XS1 - AS                ! ???
*
* Compute weights (defaults to weighted time)
*
      IF (SWEIGH.EQ.'S') THEN
         IF (TSIGFI.GT.0) THEN
            TW = 1.0/TSIGFI**2
         ELSE
            TW = TTIME * ABS(TFRES) / TTSYS**2 * 1E-6
         ENDIF
         IF (RSIGFI.GT.0) THEN
            RW = 1.0/RSIGFI**2
         ELSE
            RW = RTIME * ABS(RFRES) / RTSYS**2 * 1E-6
         ENDIF
         SW = TW + RW
*
* For SIGMA addition...
         RSIGFI = 1.0/SQRT(SW)
         TW = TW / SW
         RW = RW / SW
         SW = 1.
      ELSEIF (SWEIGH .EQ. 'T') THEN
         TW = TTIME * ABS(TFRES) / TTSYS**2 * 1E-6
         RW = RTIME * ABS(RFRES) / RTSYS**2 * 1E-6
         SW = TW + RW
         TW = TW / SW
         RW = RW / SW
         SW = 1.
      ELSE
         TW = 1.
         RW = 1.
         SW = 2.
      ENDIF
*
* Integrate
*
      RDXMAX = (AS + ABS(AR)) / 2
      RDXMIN = (AS - ABS(AR)) / 2
      TDXMAX = (AS + ABS(AT)) / 2
      TDXMIN = (AS - ABS(AT)) / 2
      IF (SALIGN.EQ.'C') GOTO 40
*
* Complex summation
      DO 30 IS = 1, SNCHAN
         XS = AS*IS + BS
         WW = 0
         SS = 0
         IF (AR.GT.0) THEN
            IRMIN = MAX(1,INT((XS-RDXMAX-BR)/AR))
            IRMAX = MIN(RNCHAN,INT((XS+RDXMAX-BR)/AR))
         ELSE
            IRMAX = MIN(RNCHAN,INT((XS-RDXMAX-BR)/AR))
            IRMIN = MAX(1,INT((XS+RDXMAX-BR)/AR))
         ENDIF
         DO 10 IR = IRMIN,IRMAX
            XR = AR*IR + BR
            DX = ABS (XR-XS)
            IF (DX .LT. RDXMAX) THEN
*
* Bad channels
               IF (RDATAY(IR).EQ.RBAD) THEN
                  IF (ALL) THEN
                     S(IS) = RBAD
                     GOTO 30
                  ENDIF
                  W = 0.
               ELSEIF (DX .LE. RDXMIN) THEN
                  W = 1.
               ELSE
                  W = 1.-(DX-RDXMIN)/(RDXMAX-RDXMIN)
               ENDIF
               WW = WW + W*RW
               SS = SS + RDATAY(IR) * W*RW
            ENDIF
10       CONTINUE
*
         IF (AT.GT.0) THEN
            ITMIN = MAX(1,INT((XS-TDXMAX-BT)/AT))
            ITMAX = MIN(TNCHAN,INT((XS+TDXMAX-BT)/AT))
         ELSE
            ITMAX = MIN(TNCHAN,INT((XS-TDXMAX-BT)/AT))
            ITMIN = MAX(1,INT((XS+TDXMAX-BT)/AT))
         ENDIF
         DO 20 IT = ITMIN, ITMAX
            XT = AT*IT + BT
            DX = ABS (XT-XS)
            IF (DX .LT. TDXMAX) THEN
               IF (TDATAY(IT).EQ.TBAD) THEN
                  IF (ALL) THEN
                     S(IS)=RBAD
                     GOTO 30
                  ENDIF
                  W = 0.
               ELSEIF (DX .LE. TDXMIN) THEN
                  W = 1.
               ELSE
                  W = 1.-(DX-TDXMIN)/(TDXMAX-TDXMIN)
               ENDIF
               WW = WW + W*TW
               SS = SS + TDATAY(IT) * W*TW
            ENDIF
20       CONTINUE
         IF (WW.NE.0) THEN
            S(IS) = SS / WW * SW
         ELSE
            S(IS) = RBAD
         ENDIF
30    CONTINUE
      GOTO 60
*
* Quick addition
40    CONTINUE
      DO 50 I = 1, SNCHAN
         IF (RDATAY(I).NE.RBAD) THEN
            S(I) = RDATAY(I)*RW
            W    = RW
         ELSEIF (ALL) THEN
            S(I) = RBAD
            GOTO 50
         ELSE
            S(I) = 0.
            W    = 0.
         ENDIF
         IF (TDATAY(I).NE.TBAD) THEN
            S(I) = S(I)+TDATAY(I)*TW
            W    = W+TW
         ELSEIF (ALL) THEN
            S(I) = RBAD
            GOTO 50
         ENDIF
         IF (W.EQ.0) THEN
            S(I) = RBAD
         ELSE
            S(I) = S(I) / W * SW
         ENDIF
50    CONTINUE
*
* Copy result in R
*
60    CONTINUE
      DO 70 I = 1, SNCHAN
70    RDATAY(I) = S(I)
      STIME = RTIME + TTIME
*
      IF (SALIGN.EQ.'V') THEN
         RVRES = AS
         RVOFF = SVOFF
         RRCHAN = (RVOFF - BS) / AS
         RFOFF = 0.
         RFRES = -RRESTF * RVRES / C
      ELSEIF (SALIGN.EQ.'F') THEN
         RFRES = AS
         RVRES = - C * RFRES / RRESTF
         RRCHAN = - BS / AS
      ENDIF
      STSYS = SQRT( STIME*ABS(RFRES)/
     $(TTIME*ABS(TFRES)/TTSYS**2 + RTIME*ABS(RFRES)/RTSYS**2) )
*
* Gestion de l'origine des spectres
      IF (RNSEQ.EQ.0) THEN
         RNSEQ = 1
         RSTART(1) = RSCAN
         REND(1) = RSCAN
      ENDIF
      IF (TNSEQ.EQ.0) THEN
         TNSEQ = 1
         TSTART(1) = TSCAN
         TEND(1) = TSCAN
      ENDIF
*
      SNSEQ = 0
      SEND(0) = -(2**14)
      IT = 1
      IR = 1
90    IF (IT.GT.TNSEQ) GOTO 91
      IF (IR.GT.RNSEQ) GOTO 92
      IF (TSTART(IT).GT.RSTART(IR)) GOTO 91
*
* TS < RS ou plus de RS
92    IF (SEND(SNSEQ).LT.TSTART(IT)-1) THEN
* nouvel intervalle
         IF(SNSEQ.EQ.MSEQ) GO TO 925
         SNSEQ = SNSEQ+1
         SSTART(SNSEQ) = TSTART(IT)
         SEND(SNSEQ) = TEND(IT)
      ELSE
* extension de l'intervalle courant
         SEND(SNSEQ) = MAX(SEND(SNSEQ),TEND(IT))
      ENDIF
      IT = IT+1
      GOTO 90
*
* RS < TS ou plus de TS
91    IF (IR.GT.RNSEQ) GOTO 93
      IF (SEND(SNSEQ).LT.RSTART(IR)-1) THEN
         IF (SNSEQ.EQ.MSEQ) GO TO 925
         SNSEQ = SNSEQ+1
         SSTART(SNSEQ) = RSTART(IR)
         SEND(SNSEQ) = REND(IR)
      ELSE
         SEND(SNSEQ) = MAX(SEND(SNSEQ),REND(IR))
      ENDIF
      IR = IR+1
      GOTO 90
*
* depassement de capacite
925   CALL MESSAGE(7,2,'ACCUMULATE','Origin table overflow')
*
* Fini
93    CONTINUE
      DO I=1,SNSEQ
         RSTART(I)=SSTART(I)
         REND(I)=SEND(I)
      ENDDO
      RNSEQ = SNSEQ
      RNUM = TNUM
      RSCAN = TSCAN
      RVER = TVER
      RTIME= STIME
      RTSYS= STSYS
      RNCHAN = SNCHAN
      LNUM = TNUM
      LVER = TVER
      LSCAN = TSCAN
      LPOSA = TAPOS
      RETURN
100   FORMAT(' R:',I5,' T:',I5)
101   FORMAT(' R:',1PG12.5,' T:',1PG12.5)
      END
*<FF>
      SUBROUTINE ADDOVER (WIN,WOUT,ERROR)
C----------------------------------------------------------------------
C SAS	Internal routine
C	Adds buffer R and T together.
C	The position compatibility checked if SMATCH is true
C	alignment with respect to frequency if SALIGN = 'F'
C	                          velocity if SALIGN = 'V'
C				  channel number if SALIGN = 'C' (default)
C
C Arguments :
C	ERROR	L	Logical error flag	Output
C----------------------------------------------------------------------
      REAL*8 C
      PARAMETER (C=299792.458)
      INCLUDE 'parameter.inc'
      INCLUDE 'index.inc'
      INCLUDE 'rdata.inc'
      INCLUDE 'structure.par'
      INCLUDE 'tdata.inc'
      INCLUDE 'structure.pat'
      INCLUDE 'setup.inc'
      INCLUDE 'work.inc'
      LOGICAL ERROR
      REAL*4 WIN (MDATA), WOUT(MDATA)
      REAL*4 S(MDATA)
      EQUIVALENCE (S,DATA_BUFFER)
      INTEGER*4 SNCHAN
*
* Variables locales
      REAL*4 DPOS2, STOL2, AT, AR, BT, BR, XT1, XT2, XR1, XR2,
     $XRMIN, XRMAX, XTMIN, XTMAX, XS1, XS2, TW, RW, SW,
     $AS, BS, RDXMAX, RDXMIN, TDXMAX, TDXMIN, XS, WW, SS,
     $XR, DX, W, XT, STIME, SVOFF, STSYS
      INTEGER IS, IR, IT, I, ITMIN, ITMAX, IRMIN, IRMAX
      LOGICAL ALL
      CHARACTER*80 CHAIN, MESS
      INTEGER*4 SNSEQ
      INTEGER*4 SSTART(0:MSEQ),SEND(0:MSEQ)
*
* Kind compatibility
*
      IF (RKIND.NE.TKIND) THEN
         CALL MESSAGE(8,3,'ADD',
     $   'Cannot add Continuum with Line data')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
* Position compatibility
*
      IF (RTYPEC.NE.TTYPEC) THEN
         CALL MESSAGE(7,3,'ACCUMULATE',
     $   'Coordinate systems are not compatible')
         ERROR = .TRUE.
         RETURN
      ELSEIF (SMATCH) THEN
         DPOS2 = ((RLAM+RLAMOF-TLAM-TLAMOF)*COS(RBET))**2
     $   + (RBET+RBETOF-TBET-TBETOF)**2
         STOL2 = STOLE ** 2
         IF (DPOS2 .GT. STOL2) THEN
            CALL MESSAGE(7,3,'ACCUMULATE',
     $      'Positions are not compatible')
            ERROR = .TRUE.
            RETURN
         ENDIF
      ENDIF
*
* Calibration Compatibility
*
      IF (SBEAMT.NE.0) THEN
         IF (ABS(RBEEFF-TBEEFF).GT.SBEAMT) THEN
            WRITE(CHAIN,101) RBEEFF,TBEEFF
            MESS = 'Different calibrations'//CHAIN
            CALL MESSAGE(8,3,'ACCUMULATE',MESS)
            ERROR = .TRUE.
            RETURN
         ENDIF
         IF (SGAINT.NE.0) THEN
            IF (ABS(RGAINI-TGAINI).GT.SGAINT) THEN
               WRITE(CHAIN,101) RBEEFF,TBEEFF
               MESS = 'Different calibrations'//CHAIN
               CALL MESSAGE(8,3,'ACCUMULATE',MESS)
               ERROR = .TRUE.
               RETURN
            ENDIF
         ENDIF
      ENDIF
*
* Frequency compatibility
*
      IF (RVTYPE.NE.TVTYPE) THEN
         MESS = 'Different velocity definition for R and T'
         CALL MESSAGE(6,2,'ACCUMULATE',MESS)
      ENDIF
      IF (SALIGN.EQ.'C') THEN
         AT = 1.
         AR = 1.
         BT = TVOFF - TVRES*TRCHAN + (RRESTF-TRESTF)/RRESTF*C
         BR = RVOFF - RVRES*RRCHAN
         IF (RRCHAN.NE.TRCHAN) THEN
            WRITE (CHAIN,101) RRCHAN,TRCHAN
            MESS = 'Different reference channels'//CHAIN
            CALL MESSAGE(6,2,'ACCUMULATE',MESS)
         ENDIF
         IF (ABS(RFRES).NE.ABS(TFRES)) THEN
            WRITE(CHAIN,101) RFRES,TFRES
            MESS = 'Different resolutions'//CHAIN
            CALL MESSAGE(6,2,'ACCUMULATE',MESS)
         ENDIF
         IF (ABS(BT-BR).GT.ABS(RVRES)*1E-2) THEN     ! SG 22-Apr-1986
            CALL MESSAGE(6,2,'ACCUMULATE',
     $      'Spectra not aligned in sky frequency')
            WRITE(6,*) BR,BT
         ENDIF
         BT = 0.
         BR = 0.
      ELSEIF ( SALIGN.EQ.'F') THEN
         AT = TFRES
         AR = RFRES
         BR = - RRCHAN * RFRES
         BT = - TRCHAN * TFRES + (TRESTF - RRESTF)
         IF (RVOFF.NE.TVOFF) THEN
            WRITE(CHAIN,101) RVOFF,TVOFF
            MESS = 'Different velocities'//CHAIN
            CALL MESSAGE(6,2,'ACCUMULATE',MESS)
         ENDIF
      ELSEIF (SALIGN.EQ.'V') THEN
         AT = TVRES
         AR = RVRES
         BT = TVOFF - TVRES*TRCHAN
         IF (RRESTF.NE.TRESTF) THEN
            WRITE(CHAIN,101) RRESTF,TRESTF
            MESS = 'Different rest frequencies'//CHAIN
            CALL MESSAGE(6,2,'ACCUMULATE',MESS)
         ENDIF
         BR = RVOFF - RVRES*RRCHAN
      ENDIF
*
* Bad channel summation
      ALL = SBAD.EQ.'O'
*
      XT1 = AT + BT
      XT2 = AT*TNCHAN + BT
      XR1 = AR + BR
      XR2 = AR*RNCHAN + BR
      XTMIN = MIN(XT1,XT2)
      XTMAX = MAX(XT1,XT2)
      XRMIN = MIN(XR1,XR2)
      XRMAX = MAX(XR1,XR2)
*
* Cas reunion seulement
      XS1 = MIN(XTMIN,XRMIN)
      XS2 = MAX(XTMAX,XRMAX)
*
* Determine parameters of sum
*
      AS = MAX(ABS(AR),ABS(AT))
      SNCHAN = (XS2 - XS1)/AS + 1.5
      IF (SNCHAN.GT.MDATA) THEN
         CALL MESSAGE(8,4,'ACCUMULATE','Too many channels requested')
         ERROR = .TRUE.
         RETURN
      ENDIF
*
      IF (AS.EQ.ABS(AR)) THEN
         XS1 = AR * NINT((XS1-BR)/AR) + BR
         SVOFF = RVOFF
      ELSE
         XS1 = AT * NINT((XS1-BT)/AT) + BT
         SVOFF = TVOFF
      ENDIF
      BS = XS1 - AS
*
* Compute weights (defaut to weighted time)
*
      IF (SWEIGH.EQ.'S') THEN
         IF (TSIGFI.GT.0) THEN
            TW = 1.0/TSIGFI**2
         ELSE
            TW = TTIME * ABS(TFRES) / TTSYS**2 * 1E-6
         ENDIF
         IF (RSIGFI.GT.0) THEN
            RW = 1.0/RSIGFI**2
         ELSE
            RW = RTIME * ABS(RFRES) / RTSYS**2 * 1E-6
         ENDIF
         SW = TW + RW
*
* For SIGMA addition...
         RSIGFI = 1.0/SQRT(SW)
*	    TW = TW / SW
*	    RW = RW / SW
*	    SW = 1.
      ELSEIF (SWEIGH .EQ. 'T') THEN
         TW = TTIME * ABS(TFRES) / TTSYS**2 * 1E-6
         RW = RTIME * ABS(RFRES) / RTSYS**2 * 1E-6
         SW = TW + RW
*	    TW = TW / SW
*	    RW = RW / SW
*	    SW = 1.
      ELSE
         TW = 1.
         RW = 1.
         SW = 2.
      ENDIF
*
* Integrate
*
      RDXMAX = (AS + ABS(AR)) / 2
      RDXMIN = (AS - ABS(AR)) / 2
      TDXMAX = (AS + ABS(AT)) / 2
      TDXMIN = (AS - ABS(AT)) / 2
      IF (SALIGN.EQ.'C') GOTO 40
*
* Complex summation
      DO 30 IS = 1, SNCHAN
         XS = AS*IS + BS
         WW = 0
         SS = 0
         IF (AR.GT.0) THEN
            IRMIN = MAX(1,INT((XS-RDXMAX-BR)/AR))
            IRMAX = MIN(RNCHAN,INT((XS+RDXMAX-BR)/AR))
         ELSE
            IRMAX = MIN(RNCHAN,INT((XS-RDXMAX-BR)/AR))
            IRMIN = MAX(1,INT((XS+RDXMAX-BR)/AR))
         ENDIF
         DO 10 IR = IRMIN,IRMAX
            XR = AR*IR + BR
            DX = ABS (XR-XS)
            IF (DX .LT. RDXMAX) THEN
*
* Bad channels
               IF (RDATAY(IR).EQ.RBAD) THEN
                  IF (ALL) THEN
                     S(IS) = RBAD
                     GOTO 30
                  ENDIF
                  W = 0.
               ELSEIF (DX .LE. RDXMIN) THEN
                  W = 1.
               ELSE
                  W = 1.-(DX-RDXMIN)/(RDXMAX-RDXMIN)
               ENDIF
               WW = WW + W*RW
               SS = SS + RDATAY(IR) * W*RW
            ENDIF
10       CONTINUE
*
         IF (AT.GT.0) THEN
            ITMIN = MAX(1,INT((XS-TDXMAX-BT)/AT))
            ITMAX = MIN(TNCHAN,INT((XS+TDXMAX-BT)/AT))
         ELSE
            ITMAX = MIN(TNCHAN,INT((XS-TDXMAX-BT)/AT))
            ITMIN = MAX(1,INT((XS+TDXMAX-BT)/AT))
         ENDIF
         DO 20 IT = ITMIN, ITMAX
            XT = AT*IT + BT
            DX = ABS (XT-XS)
            IF (DX .LT. TDXMAX) THEN
               IF (TDATAY(IT).EQ.TBAD) THEN
                  IF (ALL) THEN
                     S(IS)=RBAD
                     GOTO 30
                  ENDIF
                  W = 0.
               ELSEIF (DX .LE. TDXMIN) THEN
                  W = 1.
               ELSE
                  W = 1.-(DX-TDXMIN)/(TDXMAX-TDXMIN)
               ENDIF
               WW = WW + W*WIN(IT)
               SS = SS + TDATAY(IT) * W*WIN(IT)
            ENDIF
20       CONTINUE
         IF (WW.NE.0) THEN
            S(IS) = SS
            WOUT(IS) = WW
         ELSE
            S(IS) = RBAD
            WOUT(IS) = 0.0
         ENDIF
30    CONTINUE
      GOTO 60
*
* Quick addition
40    CONTINUE
      DO 50 I = 1, SNCHAN
         IF (RDATAY(I).NE.RBAD) THEN
            S(I) = RDATAY(I)*RW
            WOUT(I)    = RW
         ELSEIF (ALL) THEN
            S(I) = RBAD
            WOUT(I) = 0.0
            GOTO 50
         ELSE
            S(I) = 0.
            WOUT(I) = 0.
         ENDIF
         IF (TDATAY(I).NE.TBAD) THEN
            S(I) = S(I)+TDATAY(I)*WIN(I)
            WOUT(I) = WOUT(I)+WIN(I)
         ELSEIF (ALL) THEN
            S(I) = RBAD
            WOUT(I) = 0
            GOTO 50
         ENDIF
         IF (WOUT(I).EQ.0) THEN
            S(I) = RBAD
         ELSE
            S(I) = S(I)
         ENDIF
50    CONTINUE
*
* Copy result in R
*
60    CONTINUE
*
      DO I = 1, SNCHAN
         IF (WOUT(I).NE.0.0) THEN
            RDATAY(I) = S(I)/WOUT(I)
         ELSE
            RDATAY(I) = RBAD
         ENDIF
      ENDDO
      STIME = RTIME + TTIME
*
      IF (SALIGN.EQ.'V') THEN
         RVRES = AS
         RVOFF = SVOFF
         RRCHAN = (RVOFF - BS) / AS
         RFOFF = 0.
         RFRES = -RRESTF * RVRES / C
      ELSEIF (SALIGN.EQ.'F') THEN
         RFRES = AS
         RVRES = - C * RFRES / RRESTF
         RRCHAN = - BS / AS
      ENDIF
      STSYS = SQRT( STIME*ABS(RFRES)/
     $(TTIME*ABS(TFRES)/TTSYS**2 + RTIME*ABS(RFRES)/RTSYS**2) )
*
* Gestion de l'origine des spectres
      IF (RNSEQ.EQ.0) THEN
         RNSEQ = 1
         RSTART(1) = RSCAN
         REND(1) = RSCAN
      ENDIF
      IF (TNSEQ.EQ.0) THEN
         TNSEQ = 1
         TSTART(1) = TSCAN
         TEND(1) = TSCAN
      ENDIF
*
      SNSEQ = 0
      SEND(0) = -(2**14)
      IT = 1
      IR = 1
90    IF (IT.GT.TNSEQ) GOTO 91
      IF (IR.GT.RNSEQ) GOTO 92
      IF (TSTART(IT).GT.RSTART(IR)) GOTO 91
*
* TS < RS ou plus de RS
92    IF (SEND(SNSEQ).LT.TSTART(IT)-1) THEN
* nouvel intervalle
         IF(SNSEQ.EQ.MSEQ) GO TO 925
         SNSEQ = SNSEQ+1
         SSTART(SNSEQ) = TSTART(IT)
         SEND(SNSEQ) = TEND(IT)
      ELSE
* extension de l'intervalle courant
         SEND(SNSEQ) = MAX(SEND(SNSEQ),TEND(IT))
      ENDIF
      IT = IT+1
      GOTO 90
*
* RS < TS ou plus de TS
91    IF (IR.GT.RNSEQ) GOTO 93
      IF (SEND(SNSEQ).LT.RSTART(IR)-1) THEN
         IF (SNSEQ.EQ.MSEQ) GO TO 925
         SNSEQ = SNSEQ+1
         SSTART(SNSEQ) = RSTART(IR)
         SEND(SNSEQ) = REND(IR)
      ELSE
         SEND(SNSEQ) = MAX(SEND(SNSEQ),REND(IR))
      ENDIF
      IR = IR+1
      GOTO 90
*
* depassement de capacite
925   CALL MESSAGE(7,2,'ACCUMULATE','Origin table overflow')
*
* Fini
93    CONTINUE
      DO I=1,SNSEQ
         RSTART(I)=SSTART(I)
         REND(I)=SEND(I)
      ENDDO
      RNSEQ = SNSEQ
      RNUM = TNUM
      RSCAN = TSCAN
      RVER = TVER
      RTIME= STIME
      RTSYS= STSYS
      RNCHAN = SNCHAN
      LNUM = TNUM
      LVER = TVER
      LSCAN = TSCAN
      LPOSA = TAPOS
      RETURN
100   FORMAT(' R:',I5,' T:',I5)
101   FORMAT(' R:',1PG12.5,' T:',1PG12.5)
      END
*<FF>
      SUBROUTINE ADDINIT (WIN,ERROR)
C----------------------------------------------------------------------
C SAS	Internal routine
C	Initiaize COMPOSITE Adding
C	The position compatibility checked if SMATCH is true
C	alignment with respect to frequency if SALIGN = 'F'
C	                          velocity if SALIGN = 'V'
C				  channel number if SALIGN = 'C' (default)
C
C Arguments :
C	ERROR	L	Logical error flag	Output
C----------------------------------------------------------------------
      REAL*8 C
      PARAMETER (C=299792.458D0)
      INCLUDE 'rdata.inc'
      INCLUDE 'structure.par'
      INCLUDE 'setup.inc'
      INCLUDE 'work.inc'
      LOGICAL ERROR
      REAL*4 WIN (MDATA), RW
      INTEGER I
*
* Compute weights (defaut to weighted time)
*
      IF (SWEIGH.EQ.'S') THEN
         IF (RSIGFI.GT.0) THEN
            RW = 1.0/RSIGFI**2
         ELSE
            RW = RTIME * ABS(RFRES) / RTSYS**2 * 1E-6
         ENDIF
      ELSEIF (SWEIGH .EQ. 'T') THEN
         RW = RTIME * ABS(RFRES) / RTSYS**2 * 1E-6
      ELSE
         RW = 1.
      ENDIF
      DO I = 1,RNCHAN
         WIN(I) = RW
      ENDDO
      END