      PROGRAM STACKER

C     Produces time series by Fourier synthesis

      INTEGER PRTFIL, SHDFIL, FFTFIL
      PARAMETER (
     &   PRTFIL = 6, FFTFIL = 10, SHDFIL = 26,
     &   MXNSD = 500, MXNRD = 51, MXNR = 6001, MXNF = 1000,
     &   MXNT = 50000, MXN2 = 2000000,
C    &   FREQMOD = 280.0 )
     &   FREQMOD =   0.0 )

      REAL    R( MXNR ), SD( MXNSD ), RD( MXNRD ),
     &        TOUT( MXNT ), RTS( MXN2 )
      COMPLEX P( MXNR ), SNAP( MXN2 ), S( MXNF )
      CHARACTER TITLE*80, PLTTYP*10, FILNAM*80, FILEROOT*80,
     &   FILEEXT*16, PLOTTYPE*1

      EQUIVALENCE ( RTS, SNAP )

C     --- read root file name and length of name

      READ( *, * ) FILEROOT
      LAST = INDEX( FILEROOT, ' ' ) - 1

      READ( *, * ) PLOTTYPE

C     *** Read in the output times ***

      TOUT( 3 ) = -999.9
      READ( *, * ) NTOUT, ( TOUT( ITOUT ), ITOUT = 1, NTOUT )

      IF ( NTOUT .GT. MXNT ) THEN
         WRITE( PRTFIL, * ) 'NTOUT, MXNT = ', NTOUT, MXNT
         CALL ERROUT( PRTFIL, 'F', 'GETPAR',
     &      'Too many output times'  )
      ENDIF

      CALL SUBTAB( TOUT, NTOUT )
      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'Number of output times = ', NTOUT
      WRITE( PRTFIL, * ) ( TOUT( ITOUT ), ITOUT = 1, MIN( NTOUT, 51 ) )

C     *** set up source spectrum ***

      OPEN( FILE = 'FFTFIL', UNIT = FFTFIL, STATUS = 'OLD' )

      NFREQ = 0

      DO IFREQ = 1, MXNF
         READ( FFTFIL, *, END = 1000 ) F, SR, SI
         WRITE( *, * ) F, SR, SI
         S( IFREQ ) = CMPLX( SR, SI )
         NFREQ = IFREQ
      END DO

 1000 CONTINUE

      FILNAM = 'work/' // FILEROOT(1:LAST) // '1.shd'

      CALL OPENFILES( FILNAM, PLOTTYPE, SD, NSD, MXNSD,
     &    RD, NRD, MXNRD, R, NR, MXNR, TOUT, NTOUT )

      NSTORE = NTOUT * NRD * NR
      IF ( NSTORE .GT. MXN2 ) THEN
         WRITE( PRTFIL, * ) 'NSTORE, MXN2 = ', NSTORE, MXN2
         CALL ERROUT( PRTFIL, 'F', 'GETPAR',
     &      'Too much field data'  )
      ENDIF

C     *** Main loop ***

      DO ISD = 1, NSD

C        --- zero out the matrices
         DO I = 1, NTOUT * NRD * NR
            SNAP( I ) = 0.0
         END DO

         DO I = 1, NTOUT * NRD
            RTS( I ) = 0.0
         END DO

C        *** frequency loop ***

         DO IFREQ = 1, NFREQ

            IF      ( IFREQ .LE. 9   ) THEN
               WRITE( FILEEXT, 2000 ) IFREQ
 2000          FORMAT( I1.1, '.shd' )
            ELSE IF ( IFREQ .LE. 99  ) THEN
               WRITE( FILEEXT, 3000 ) IFREQ
 3000          FORMAT( I2.2, '.shd' )
            ELSE IF ( IFREQ .LE. 999 ) THEN
               WRITE( FILEEXT, 4000 ) IFREQ
 4000          FORMAT( I3.3, '.shd' )
            ELSE
               WRITE( FILEEXT, 5000 ) IFREQ
 5000          FORMAT( I4.4, '.shd' )
            END IF

            FILNAM = 'work/' // FILEROOT(1:LAST) // FILEEXT
            WRITE( *, * ) FILNAM(1:40)

C           *** Read header ***

            CALL RDHEAD( SHDFIL, FILNAM, MXNR, MXNSD, MXNRD,
     &         FREQ, R, NR, DELTAR, PLTTYP, XS, YS, THETA,
     &         SD, NSD, RD, NRD, TITLE )

            FREQ = FREQ - FREQMOD   ! frequency downshift (demodulation)

C           *** read in a single frequency frame ***

            DO IRD = 1, NRD
               IREC =  6 + NRD * ( ISD - 1 ) + IRD

               READ( SHDFIL, REC = IREC ) ( P( IR ), IR = 1, NR )

C              --- add spectral component to matrix 

               CALL ADDIN( IRD, FREQ, S( IFREQ ), P, RTS, SNAP,
     &                   NRD, NR, TOUT, NTOUT, PLOTTYPE )

            END DO   ! next receiver depth IRD

            CLOSE( SHDFIL )

         END DO   ! next frequency IFREQ

C        *** write results ***

         CALL WRTDAT( RTS, SNAP, TOUT, RD, NTOUT, NRD, NR, PLOTTYPE )

      END DO   ! next ISD

      STOP
      END

C**********************************************************************C

      SUBROUTINE OPENFILES( FILNAM, PLOTTYPE, SD, NSD, MXNSD,
     &   RD, NRD, MXNRD, R, NR, MXNR, TOUT, NTOUT )

C     Opens up the output files and writes header data

      INTEGER SHDFIL, RTSFIL
      PARAMETER ( SHDFIL = 26, RTSFIL = 35 )

      REAL      R( * ), SD( * ), RD( * ), TOUT( * )
      CHARACTER TITLE*80, PLTTYP*10, FILNAM*80, PLOTTYPE*1

C     *** Read in header info from first frequency file ***

      CALL RDHEAD( SHDFIL, FILNAM, MXNR, MXNSD, MXNRD,
     &      FREQ, R, NR, DELTAR, PLTTYP, XS, YS, THETA,
     &      SD, NSD, RD, NRD, TITLE )

      CLOSE( SHDFIL )

C    *** Open output file ***

      IF ( PLOTTYPE .EQ. 'D' ) THEN   ! *** depth stack ***
         OPEN ( FILE = 'RTSFIL', UNIT = RTSFIL, STATUS = 'NEW',
     &          FORM = 'FORMATTED' )

         WRITE( RTSFIL, * ) '''' // TITLE(1:75) // ''''
         WRITE( RTSFIL, * ) NRD, ( RD( IR ), IR = 1, NRD )

      ELSE IF ( PLOTTYPE .EQ. 'R' ) THEN   ! *** range stack ***
         OPEN ( FILE = 'RTSFIL', UNIT = RTSFIL, STATUS = 'NEW',
     &          FORM = 'FORMATTED' )

         WRITE( RTSFIL, * ) '''' // TITLE(1:75) // ''''
         WRITE( RTSFIL, * ) NR, ( R( IR ), IR = 1, NR )

      ELSE   ! *** snapshot ***

         CALL HEADER( TITLE, TOUT, NTOUT, RD, NRD,
     &      R, NR, FREQ, ATTEN, PLTTYP, XS, YS, THETA )

      END IF

      RETURN
      END

C**********************************************************************C

      SUBROUTINE ADDIN( IRD, FREQ, S, P, RTS, SNAP, NRD, NR,
     &    TOUT, NTOUT, PLOTTYPE )

      COMPLEX CI
      PARAMETER ( CI = ( 0.0, 1.0 ), TWOPI = 2.0 * 3.141592 )

      REAL    TOUT( * ), RTS( NTOUT, * )
      COMPLEX P( * ), SNAP( NR, NTOUT, NRD ), S, W
      CHARACTER PLOTTYPE*1

C     Add in the spectral component
C     Note optimal loop order is debatable
C     Since RTS is real, we are taking the real part of the transform ...

      IF ( PLOTTYPE .EQ. 'D' ) THEN   ! * depth stack *
         DO ITOUT = 1, NTOUT
             RTS( ITOUT, IRD ) = RTS( ITOUT, IRD ) +
     &          S * P( 1 ) * EXP( CI * TWOPI * FREQ * TOUT( ITOUT ) )
         END DO

      ELSE IF ( PLOTTYPE .EQ. 'R') THEN   ! * range stack *
         DO ITOUT = 1, NTOUT
            W = S * EXP( CI * TWOPI * FREQ * TOUT( ITOUT ) )
            DO IR = 1, NR
               RTS( ITOUT, IR ) = RTS( ITOUT, IR ) + W * P( IR )
            END DO
         END DO

      ELSE   ! * snapshot *
         DO ITOUT = 1, NTOUT
            W = S * EXP( CI * TWOPI * FREQ * TOUT( ITOUT ) )
            DO IR = 1, NR
               SNAP( IR, ITOUT, IRD ) = SNAP( IR, ITOUT, IRD ) +
     &           W * P( IR )
            END DO
         END DO
      END IF

      RETURN
      END

C**********************************************************************C

      SUBROUTINE WRTDAT( RTS, SNAP, TOUT, RD, NTOUT, NRD, NR, PLOTTYPE )

C     Writes data to disk

      INTEGER SNAFIL, RTSFIL
      PARAMETER ( SNAFIL = 25, RTSFIL = 35 )

      REAL      TOUT( * ), RTS( NTOUT, * ), RD( * )
      COMPLEX   SNAP( NR, NTOUT, NRD )
      CHARACTER PLOTTYPE*1

      IF ( PLOTTYPE .EQ. 'D' ) THEN   ! depth stack
         DO ITOUT = 1, NTOUT
            WRITE( RTSFIL, * ) TOUT( ITOUT ),
     &         ( RTS( ITOUT, IRD ), IRD = 1, NRD )
         END DO

      ELSE IF ( PLOTTYPE .EQ. 'R') THEN   ! range stack
         DO ITOUT = 1, NTOUT
            WRITE( RTSFIL, * ) TOUT( ITOUT ),
     &         ( RTS( ITOUT, IR ), IR = 1, NR )
         END DO

      ELSE   ! snapshot
         DO ITOUT = 1, NTOUT
            DO IRD = 1, NRD
               IREC =  6 + NRD * ( ITOUT - 1 ) + IRD
               WRITE( SNAFIL, REC = IREC )
     &                 ( SNAP( IR, ITOUT, IRD ), IR = 1, NR )
            END DO
         END DO
      END IF

      RETURN
      END

