      PROGRAM SPARC

C     Finite element time domain fast-field program
C     Michael B. Porter   10/87

      INCLUDE 'COMSPA.FOR'

      REAL*8 TEND
      COMPLEX GREEN( 300000 )
      REAL RTSRD( 300000 ), RTSRR( 300000 )
      EQUIVALENCE ( GREEN, RTSRD, RTSRR )

      CALL CLTIME

      CALL GETPAR

      DO MED = 1, NMEDIA
         H( MED ) = ( DEPTH( MED + 1 ) - DEPTH( MED ) ) / N( MED )
      END DO

      CALL INIT
      CALL WRTHED
      CALL KERNEL( GREEN, RTSRD, RTSRR )

      CALL RDTIME( TEND )
      WRITE( PRTFIL, 5000) TEND
 5000 FORMAT(' CPU TIME: ', G15.5)

      STOP
      END
C**********************************************************************C
      SUBROUTINE GETPAR

C     Read in the ENVFIL data

      INCLUDE 'COMSPA.FOR'

      REAL*8 FREQ, RHOTD, RHOBD, BUMDEN, ETA, XI
      COMPLEX*16 CPTD, CSTD, CPBD, CSBD

C     *** Call READIN to get most of the environmental data ***

      TITLE = 'SPARC-   '
      CALL READIN( TITLE, FREQ, MXMED, NMEDIA,
     &   TOPOPT, CPTD, CSTD, RHOTD, BUMDEN, ETA, XI, N, SIGMA, DEPTH, 
     &   BOTOPT, CPBD, CSBD, RHOBD,
     &   ENVFIL, PRTFIL )

C     *** Spectral limits ***

      READ(  ENVFIL, *    ) CLOW, CHIGH
      WRITE( PRTFIL, 6000 ) CLOW, CHIGH
 6000 FORMAT( /, ' CLOW = ', G12.5, '  CHIGH = ', G12.5 )

C     *** Maximum range for calculations ***

      READ(  ENVFIL, * ) RMAX
      WRITE( PRTFIL, * ) 'RMAX = ', RMAX

C     *** Read source/receiver depths ***

      ZMIN = DEPTH( 1 )
      ZMAX = DEPTH( NMEDIA + 1 )

      CALL SDRD( ENVFIL, PRTFIL, ZMIN, ZMAX,
     &   MAXNSD, NSD, SD,
     &   MAXNRD, NRD, RD )

      OMEGA2 = ( 2.0 * PI * FREQ ) ** 2

      CPT  = CPTD
      CST  = CSTD
      RHOT = RHOTD
      CPB  = CPBD
      CSB  = CSBD
      RHOB = RHOBD

C     *  Upper and lower frequency limits *

      READ( ENVFIL, * ) PULSE, FMIN, FMAX

C     ------ integration parameters
      RKMIN = 2.0 * PI * FMIN / CHIGH
      RKMAX = 2.0 * PI * FMAX / CLOW
      IF ( CHIGH .GE. 1.0E6 ) RKMIN = 0.0

      NKPTS = 1000.0 * RMAX * ( RKMAX - RKMIN ) / ( 2.0 * PI )

      IF ( NKPTS .GT. MXNK ) THEN
         WRITE( PRTFIL, * ) 'NKPTS, MXNK = ', NKPTS, MXNK
         CALL ERROUT( PRTFIL, 'F', 'GETPAR',
     &      'Too many pts in k-space'  )
      ENDIF

      WRITE( PRTFIL, * ) 'NKPTS = ', NKPTS

      DELTAK = ( RKMAX - RKMIN ) / ( NKPTS - 1 )
      DO IK = 1, NKPTS
         RK( IK ) = RKMIN + ( IK - 1 ) * DELTAK
      END DO

      IF ( PULSE(1:1) .EQ. 'P' ) THEN
         WRITE( PRTFIL, * ) 'Pseudo-gaussian pulse'
      ELSE IF ( PULSE(1:1) .EQ. 'R' ) THEN
         WRITE( PRTFIL, * ) 'Ricker wavelet'
      ELSE IF ( PULSE(1:1) .EQ. 'A' ) THEN
         WRITE( PRTFIL, * ) 'Approximate Ricker wavelet'
      ELSE IF ( PULSE(1:1) .EQ. 'S' ) THEN
         WRITE( PRTFIL, * ) 'Single sine source'
      ELSE IF ( PULSE(1:1) .EQ. 'H' ) THEN
         WRITE( PRTFIL, * ) 'Hanning weighted four sine pulse'
      ELSE IF ( PULSE(1:1) .EQ. 'N' ) THEN
         WRITE( PRTFIL, * ) 'N-wave pulse'
      ELSE IF ( PULSE(1:1) .EQ. 'M' ) THEN
         WRITE( PRTFIL, * ) 'Miracle-wave pulse'
      ELSE IF ( PULSE(1:1) .EQ. 'G' ) THEN
         WRITE( PRTFIL, * ) 'Gaussian pulse'
      ELSE IF ( PULSE(1:1) .EQ. 'F' ) THEN
         WRITE( PRTFIL, * ) 'Source time series from File'
      ELSE IF ( PULSE(1:1) .EQ. 'B' ) THEN
         WRITE( PRTFIL, * ) 'Source time series reversed from file'
      ELSE
         CALL ERROUT( PRTFIL, 'F', 'GETPAR',
     &      'Unknown source type'  )
      ENDIF

C     *** Read receiver ranges ***

      RR( 3 ) = -999.9
      READ( ENVFIL, * ) NRR, ( RR( IRR ), IRR = 1, NRR )

      IF ( NRR .GT. MXNRR ) THEN
         WRITE( PRTFIL, * ) 'NRR, MXNRR = ', NRR, MXNRR
         CALL ERROUT( PRTFIL, 'F', 'GETPAR',
     &      'Too many receiver ranges'  )
      ENDIF

      CALL SUBTAB( RR, NRR )
      WRITE( PRTFIL, * )
      WRITE( PRTFIL, * ) 'Number of receiver ranges = ', NRR
      WRITE( PRTFIL, * ) ( RR( IRR ), IRR = 1, MIN( NRR, 51 ) )

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

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

      IF ( NTOUT .GT. MXNTOUT ) THEN
         WRITE( PRTFIL, * ) 'NTOUT, MXNTOUT = ', NTOUT, MXNTOUT
         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     *** Integration parameters ***
C     ALPHA = 0.0  lumped mass matrix,
C             1.0  consistent mass matrix
C     BETA = 0.0  standard explicit,
C            0.25 standard implicit 
 
      READ(  ENVFIL, * ) TSTART, TMULT, ALPHA, BETA, V
      WRITE( PRTFIL, * ) 'TSTART = ', TSTART
      WRITE( PRTFIL, * ) 'TMULT  = ', TMULT
      WRITE( PRTFIL, * ) 'ALPHA  = ', ALPHA
      WRITE( PRTFIL, * ) 'BETA   = ', BETA
      WRITE( PRTFIL, * ) 'V      = ', V

      CLOSE ( ENVFIL )
      RETURN
      END
C**********************************************************************C
      SUBROUTINE INIT

C     Initializes arrays defining difference equations

      INCLUDE 'COMSPA.FOR'

      REAL*8 FREQ
      COMPLEX*16 CP( MXN ), CS( MXN )

      CMIN   =  1.0E6
      CMAX   = -1.0E6
      CROSST =  1.0E6
      LOC( 1 ) = 0
      J = 0

C     *** Loop over successive media ***

      DO MED = 1, NMEDIA

         IF ( SIGMA( MED ) .NE. 0.0 ) THEN
            CALL ERROUT( PRTFIL, 'F', 'INIT',
     &         'Rough interfaces not allowed'  )
         ENDIF

         IF ( MED .NE. 1 ) LOC(MED) = LOC(MED-1) + N(MED-1) + 1
         N1 = N( MED ) + 1

         CALL PROFIL( DEPTH, CP, CS, RHO( LOC( MED ) + 1 ),
     &      MED, N1, FREQ, TOPOPT(1:1), TOPOPT(3:3), 'TAB',
     &      ENVFIL, PRTFIL )

         DO I = 1, N1
            CPR  = REAL( CP( I ) )
            CMIN = MIN( CPR, CMIN )
            CMAX = MAX( CPR, CMAX )

            CROSST1 = REAL( H( MED ) / CP( I ) )
            CROSST  = AMIN1( CROSST1, CROSST )

            J = J + 1
            CPT = CP( I )   ! convert to single precision
            C2R( J ) =  REAL( CPT**2 )  ! CHECK THIS
            C2I( J ) = AIMAG( CPT**2 ) / SQRT( OMEGA2 ) / C2R( J )
         END DO

      END DO

C     *** Tabulate z-coordinates ***

      J = 1
      Z( 1 ) = DEPTH( 1 )
      NTOT = 0

      DO MED = 1, NMEDIA
         DO I = 1, N( MED )
            J = J + 1
            Z( J ) = Z( J - 1 ) + H( MED )
         END DO
         NTOT = NTOT + N( MED )
      END DO

      NTOT1 = NTOT + 1

C     * Compute weights for source/rcvr depth interpolation *

      CALL WEIGHT( Z, NTOT1, SD, NSD, WS, ISD )
      CALL WEIGHT( Z, NTOT1, RD, NRD, WR, IRD )

      RETURN
      END
C**********************************************************************C
      SUBROUTINE KERNEL( GREEN, RTSRD, RTSRR )

C     Solve system for a sequence of k-values

      INCLUDE 'COMSPA.FOR'

      REAL RTSRD( NRD, * ), RTSRR( NRR, * )
      COMPLEX GREEN( NTOUT, NRD, * )

C     ------ Clear RTS

      IF ( TOPOPT(4:4) .EQ. 'D' ) THEN
         DO ITOUT = 1, NTOUT
            DO IR = 1, NRD
               RTSRD( IR, ITOUT ) = 0.0
            END DO
         END DO
      END IF

      IF ( TOPOPT(4:4) .EQ. 'R' ) THEN
         DO ITOUT = 1, NTOUT
            DO IR = 1, NRR
               RTSRR( IR, ITOUT ) = 0.0
            END DO
         END DO
      END IF

C     * Courant condition to set time step *

      DELT = TMULT /
     &       SQRT( 1.0 / CROSST**2 + ( 0.5 * CMAX * RK(NKPTS) ) **2 )

      WRITE( PRTFIL, * ) 'Time step = ', DELT
      WRITE( PRTFIL, * ) 'Estimated fl. pt. ops (millions) = ',
     &   ( TOUT( NTOUT ) / DELT ) * NKPTS * NTOT1 / 25000

C     *** Loop over spectral components ***

      IKCORE = 0
      IFLUSH = 0      ! Keeps track of number of records flushed

      DO IK = 1, NKPTS
         IKCORE = IKCORE + 1
         X = RK( IK ) ** 2
C         DELT = TMULT /
C     &             SQRT( 1.0 / CROSST**2 + ( 0.5 * CMAX * X ) **2 )
C         IF ( 10*(IK/10) .EQ. IK ) WRITE( 6, * ) 'IK, NKPTS', IK, NKPTS
         WRITE( 6, * ) 'IK, NKPTS', IK, NKPTS

C        *** March that component for all time ***

         CALL MARCH( X, GREEN, IKCORE, RTSRD, RTSRR )

C        *** Snapshot: need to flush buffer?  ***

         IF ( TOPOPT(4:4) .EQ. 'S' ) THEN
         IF ( IKCORE .EQ. NKCORE .OR. IK .EQ. NKPTS ) THEN
C           ------ number of records per S/R combo:
            NRECSR = ( NKPTS - 1 ) / NKCORE + 1

            DO ITOUT = 1, NTOUT
               DO IR = 1, NRD
C                 ------ Index of source/rcvr combo
                  IG = ( ITOUT - 1 ) * NRD + IR
                  IREC = 6 + IG * NRECSR + IFLUSH
                  WRITE( GRNFIL, REC = IREC)
     &               ( GREEN( ITOUT, IR, IKCORE ), IKCORE = 1, NKCORE )
               END DO
            END DO

            IKCORE = 0
            IFLUSH = IFLUSH + 1
         ENDIF
         ENDIF
      END DO

C     *** Write out RTS (vertical array) *

      IF ( TOPOPT(4:4) .EQ. 'D' ) THEN
         SCALE = 1.0 / SQRT( PI * 1000.0 * RR( 1 ) )
         DO ITOUT = 1, NTOUT
            WRITE( RTSFIL, * ) TOUT( ITOUT ),
     &         ( SCALE * RTSRD( IR, ITOUT ), IR = 1, NRD )
 5500       FORMAT( 12G15.6 )
         END DO
      ENDIF

C     *** Write out RTS (horizontal array) ***

      IF ( TOPOPT(4:4) .EQ. 'R' ) THEN
         DO ITOUT = 1, NTOUT
            WRITE( RTSFIL, * ) TOUT( ITOUT ),
     &         ( RTSRR( IRR, ITOUT ), IRR = 1, NRR )
         END DO
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE WRTHED

C     Writes headers for disk files

      INCLUDE 'COMSPA.FOR'
      REAL*8 FREQ

      IF ( TOPOPT(4:4) .EQ. 'S' ) THEN

C        *** Case of a snapshot ***

         NSR    = NTOUT * NRD
         NKCORE = MIN( 300000 / NSR, NKPTS, 4094 )
         FREQ   = 0.0
         CALL HEADER( FREQ )

      ELSE IF ( TOPOPT(4:4) .EQ. 'R' ) THEN

C        *** Horizontal array ***

         OPEN ( FILE = 'RTSFIL', UNIT = RTSFIL, STATUS = 'NEW',
     &          FORM = 'FORMATTED' )

         WRITE( RTSFIL, * ) '''' // TITLE(1:75) // ''''
         WRITE( RTSFIL, * ) NRR, ( RR( IRR ), IRR = 1, NRR )
 5500    FORMAT( I5, 11G15.6 )

      ELSE IF ( TOPOPT(4:4) .EQ. 'D' ) THEN

C        *** Vertical array ***

         OPEN ( FILE = 'RTSFIL', UNIT = RTSFIL, STATUS = 'NEW',
     &          FORM = 'FORMATTED' )

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

      RETURN
      END
C**********************************************************************C
      SUBROUTINE HEADER( FREQ )

C     Routine to write header to disk file

      INCLUDE 'COMSPA.FOR'

      REAL*8 FREQ
      CHARACTER PLTTYP*10

      LRECL  = MAX( 40, 2 * NKCORE, NTOUT + 1, NRD + 1 )
      NRECSR = ( NKPTS - 1 ) / NKCORE + 1
      ATTEN  = 0.0
      PLTTYP = 'GREEN'
      XS = 0.0
      YS = 0.0
      THETA = 0.0

      OPEN ( FILE = 'GRNFIL', UNIT = GRNFIL, STATUS = 'NEW',
     &   ACCESS = 'DIRECT', RECL = 4 * LRECL, FORM = 'UNFORMATTED')

      WRITE( GRNFIL, REC = 1 ) LRECL, TITLE
      WRITE( GRNFIL, REC = 2 ) PLTTYP, XS, YS, THETA
      WRITE( GRNFIL, REC = 3 ) REAL( FREQ ), NTOUT, NRD, NKPTS,
     &   NKCORE, DELTAK, ATTEN, NRECSR
      WRITE( GRNFIL, REC = 4 ) ( TOUT( I ), I = 1, NTOUT )
      WRITE( GRNFIL, REC = 5 ) ( RD(   I ), I = 1, NRD )

      DO IREC = 1, NRECSR
         ISTART = ( IREC - 1 ) * NKCORE
         WRITE( GRNFIL, REC = 5 + IREC)
     &      ( RK( ISTART + I ), I = 1, NKCORE)
      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE MARCH( X, GREEN, IKCORE, RTSRD, RTSRR )

C     March a single spectral component forward in time

      INCLUDE 'COMSPA.FOR'

      REAL RTSRD( NRD, * ), RTSRR( NRR, * )
      COMPLEX AD2( MXN ), AE2( MXN ), AD1( MXN ), AE1( MXN ),
     &        AD0( MXN ), AE0( MXN ),
     &   RV1( MXN ), RV2( MXN ), RV3( MXN ), RV4( MXN ),
     &   U0( MXN ), U1( MXN ), U2( MXN ), GREEN( NTOUT, NRD, * )

C     *** Assemble mass and stiffness matrices ***

      RKT  = SQRT( X )
      NODE = 1
      L    = LOC( 1 ) + 1

      AD2( 1 ) = 0.0
      AD1( 1 ) = 0.0
      AD0( 1 ) = 0.0

      DO MED = 1, NMEDIA
         HELT = H( MED )

         DO I = 1, N(MED)
            RHOELT = ( RHO( L ) + RHO( L + 1 ) ) / 2.0
            C2RELT = ( C2R( L ) + C2R( L + 1 ) ) / 2.0
            C2IELT = ( C2I( L ) + C2I( L + 1 ) ) / 2.0

            CALL CONTRIB( NODE, X, RKT, V, RHOELT, C2RELT, C2IELT,
     &         HELT, DELT, AD2, AE2, AD1, AE1, AD0, AE0, ALPHA, BETA )

            NODE = NODE + 1
            L    = L    + 1
         END DO

         L = L + 1
      END DO

C     * Factor A2 *

      CALL FACTOR( NTOT1, AD2, AE2, IERR, RV1, RV2, RV3, RV4 )

C     * Initialize pressure vectors *

      DO NODE = 1, NTOT1
         U0( NODE ) = 0.0
         U1( NODE ) = 0.0
         U2( NODE ) = 0.0
      END DO

C     *** Some initialization ***

      IF ( PULSE(4:4) .EQ. 'L' .OR. PULSE(4:4) .EQ. 'B' ) THEN
         FLOCUT = RKT * CLOW / ( 2.0 * PI )
      ELSE
         FLOCUT = 0.0
      ENDIF

      IF ( PULSE(4:4) .EQ. 'H' .OR. PULSE(4:4) .EQ. 'B' ) THEN
         FHICUT = RKT * CHIGH / ( 2.0 * PI )
      ELSE
         FHICUT = 10.0 * FMAX
      ENDIF

      TIME   = 0.0

C     *** Begin forward march ***

      ITOUT = 1
      DO ITIME = 1, MXIT
         TIME = TSTART + ( ITIME - 1 ) * DELT

C        * Take a step forward in time *

         CALL STEP( AD0, AE0, AD1, AE1, AD2, AE2, RV1, RV2, RV3, RV4,
     &      U0, U1, U2, TIME, RKT, FLOCUT, FHICUT )

C        * Extract soln for desired receivers *

         CALL EXTRACT( U0, U1, GREEN, IKCORE, RTSRD, RTSRR,
     &      TIME, ITOUT, RKT )

         IF ( ITOUT .GT. NTOUT ) RETURN
      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE CONTRIB( NODE, X, RKT, V, RHOELT, C2RELT, C2IELT,
     &   HELT, DELT, AD2, AE2, AD1, AE1, AD0, AE0, ALPHA, BETA )

C     Computes the contribution for the given element

      COMPLEX AD2( * ), AE2( * ), AD1( * ), AE1( * ), AD0( * ),
     &   AE0( * ), KE( 2, 2 ), CE( 2, 2 ), ME( 2, 2 ), Y, Z, CI
      PARAMETER ( CI = ( 0.0, 1.0 ) )

C     ALPHA controls lumping and BETA controls implicitness

      DELT2 = DELT ** 2
      RHOH  = RHOELT * HELT

C     * Elemental stiffness matrix ( +k2 term ) *

      Y = 1.0 + CI * RKT * V * C2IELT
      Z = X * ( Y - V * V / C2RELT )
      KE( 1, 1 ) =  Y / RHOH + (3.0-ALPHA) * HELT * Z / RHOELT / 6.0
      KE( 1, 2 ) = -Y / RHOH +      ALPHA  * HELT * Z / RHOELT / 6.0
      KE( 2, 2 ) =  Y / RHOH + (3.0-ALPHA) * HELT * Z / RHOELT / 6.0

C     * Elemental damping matrix *

      Z = 2.0*CI*RKT*V/C2RELT
      CE( 1, 1 ) = C2IELT * ( 1.0/RHOH + (3.0-ALPHA)*HELT*X/RHOELT/6.0)
     &                                 + (3.0-ALPHA)*HELT*Z/RHOELT/6.0
      CE( 1, 2 ) = C2IELT * (-1.0/RHOH +      ALPHA *HELT*X/RHOELT/6.0)
     &                                 +      ALPHA *HELT*Z/RHOELT/6.0
      CE( 2, 2 ) = C2IELT * ( 1.0/RHOH + (3.0-ALPHA)*HELT*X/RHOELT/6.0)
     &                                 + (3.0-ALPHA)*HELT*Z/RHOELT/6.0

C     * Elemental mass matrix *

      ME( 1, 1 ) = (3.0-ALPHA) * HELT / ( RHOELT * C2RELT ) / 6.0
      ME( 1, 2 ) =      ALPHA  * HELT / ( RHOELT * C2RELT ) / 6.0
      ME( 2, 2 ) = (3.0-ALPHA) * HELT / ( RHOELT * C2RELT ) / 6.0

C     * A2 matrix *

      AD2( NODE   ) = AD2( NODE ) +
     &                ME(1,1) + 0.5*DELT*CE(1,1) + BETA*DELT2*KE(1,1)
      AE2( NODE+1 ) = ME(1,2) + 0.5*DELT*CE(1,2) + BETA*DELT2*KE(1,2)
      AD2( NODE+1 ) = ME(2,2) + 0.5*DELT*CE(2,2) + BETA*DELT2*KE(2,2)

C     * A1 matrix *

      AD1( NODE   ) = AD1( NODE ) +
     &                2.0*ME(1,1) - ( 1.0 - 2.0*BETA )*DELT2*KE(1,1)
      AE1( NODE+1 ) = 2.0*ME(1,2) - ( 1.0 - 2.0*BETA )*DELT2*KE(1,2)
      AD1( NODE+1 ) = 2.0*ME(2,2) - ( 1.0 - 2.0*BETA )*DELT2*KE(2,2)

C     * A0 matrix *

      AD0( NODE   ) = AD0( NODE )
     &                -ME(1,1) + 0.5*DELT*CE(1,1) - BETA*DELT2*KE(1,1)
      AE0( NODE+1 ) = -ME(1,2) + 0.5*DELT*CE(1,2) - BETA*DELT2*KE(1,2)
      AD0( NODE+1 ) = -ME(2,2) + 0.5*DELT*CE(2,2) - BETA*DELT2*KE(2,2)

      RETURN
      END
C**********************************************************************C
      SUBROUTINE STEP( AD0, AE0, AD1, AE1, AD2, AE2,
     &   RV1, RV2, RV3, RV4, U0, U1, U2, TIME, RKT, FLOCUT, FHICUT )

C     Take a time step

C     This is the inner loop.  Note that a significant speed-up can
C     be obtained by using real arithmetic and eliminating the
C     tridiagonal solver if the Hilbert transform is bypassed and the
C     explicit solver is used.

      INCLUDE 'COMSPA.FOR'

      COMPLEX AD2( * ), AE2( * ), AD1( * ), AE1( * ), AD0( * ),
     &   AE0( * ),      RV1( * ), RV2( * ), RV3( * ), RV4( * ),
     &   U0( * ), U1( * ), U2( * ), ST( MAXNSD )
      CHARACTER PULSTITL*60

      OMEGA = SQRT( OMEGA2 )
      DELT2 = DELT ** 2

C     *** Form U2TEMP = A1*U1 + A0*U0 + S ***

C     ------ Surface point
      J = 1
      U2(J) = AD1(J) * U1(J) + AE1(J+1) * U1(J+1) +
     &        AD0(J) * U0(J) + AE0(J+1) * U0(J+1)

C     ------ Interior points
      DO J = 2, NTOT1 - 1
         U2(J) = AD1(J) * U1(J) + AE1(J) * U1(J-1) + AE1(J+1) * U1(J+1)
     &         + AD0(J) * U0(J) + AE0(J) * U0(J-1) + AE0(J+1) * U0(J+1)
      END DO

C     ------ Bottom point
      J = NTOT1
      U2(J) = AD1(J) * U1(J) + AE1(J) * U1(J-1) + 
     &        AD0(J) * U0(J) + AE0(J) * U0(J-1)

C     *** Source terms ***

      NTPTS = 1
      CALL SOURCE( TIME, ST, SD, NSD, NTPTS,
     &             OMEGA, FLOCUT, FHICUT, PULSE, PULSTITL )

      MED = 1	  ! Assumes source in first medium

      DO IS = 1, NSD
         ISDT = ISD( IS )
         JS   = LOC( 1 ) + ISDT + ( MED - 1 )

         ST( IS ) = DELT2 * ST( IS ) * EXP( -CI * RKT * V * TIME )

         U2( JS   ) = U2( JS   ) + ( 1.0 - WS( IS ) ) * ST( IS )
         U2( JS+1 ) = U2( JS+1 ) +         WS( IS )   * ST( IS )
      END DO

C     *** Solve A2*U2 = U2TEMP ***

C     --- Implicit or explicit solver?

      IF ( ALPHA .EQ. 0.0 .AND. BETA .EQ. 0.0 ) THEN
         DO I = 1, NTOT1
            U2( I ) = U2( I ) / AD2( I )
         END DO
      ELSE
         CALL BACKSB( NTOT1, AD2, AE2, IERR, RV1, RV2, RV3, RV4, U2 )
      ENDIF

C     Do a roll

      DO J = 1, NTOT1
         U0( J ) = U1( J )
         U1( J ) = U2( J )
      END DO

C     * Boundary conditions (natural is natural) *

      IF ( TOPOPT(2:2) .EQ. 'V' ) U1( 1 )     = 0.0
      IF ( BOTOPT(1:1) .EQ. 'V' ) U1( NTOT1 ) = 0.0

      RETURN
      END
C**********************************************************************C
      SUBROUTINE EXTRACT( U0, U1, GREEN, IKCORE, RTSRD, RTSRR,
     &   TIME, ITOUT, RKT )

C     Extract solution (snapshot, vertical or horizontal array)

      INCLUDE 'COMSPA.FOR'

      REAL RTSRD( NRD, * ), RTSRR( NRR, * )
      COMPLEX U0( * ), U1( * ), UT1, UT2, U, CONST, REST,
     &   GREEN( NTOUT, NRD, * )

C     *** Three cases: snapshot, vertical, horizontal time series ***

C     simulate:
C     DO WHILE ( TIME + DELT .GE. TOUT( ITOUT ) .AND. ITOUT .LE. NTOUT )

      DO IDUM = 1, 1000
         IF ( TIME + DELT .LT. TOUT( ITOUT ) .OR.
     &      ITOUT .GT. NTOUT ) RETURN

C        ------ Weight for temporal interpolation:
         WT = ( TOUT( ITOUT ) - TIME ) / DELT

         IF ( TOPOPT(4:4) .EQ. 'S' ) THEN

C           *** Case of a snapshot ***

            CONST = EXP( CI * RKT * V * TOUT( ITOUT ) )

            DO IR = 1, NRD
C              ------ Linear interpolation in depth
               I = IRD( IR )
               UT1 = U0( I ) + WR( IR ) * ( U0( I+1 ) - U0( I ) )
               UT2 = U1( I ) + WR( IR ) * ( U1( I+1 ) - U1( I ) )
C              ------ Linear interpolation in time
               GREEN( ITOUT, IR, IKCORE ) = CONST*(UT1 + WT*(UT2-UT1))
            END DO

         ELSE IF ( TOPOPT(4:4) .EQ. 'D' ) THEN

C           *** Case of an RTS (vertical array) ***

            RM = 1000.0 * RR( 1 )     ! Range (in meters)
C           CONST = -SQRT(2.0) * DELTAK * SQRT(RKT) * COS(RKT*RM-PI/4.0)
            CONST = SQRT( 2.0 ) * DELTAK * SQRT( RKT ) *
     &         EXP( CI*( RKT * ( V * TOUT( ITOUT ) - RM ) + PI / 4.0 ) )

            DO IR = 1, NRD
C              ------ Linear interpolation in depth
               I = IRD( IR )

               UT1 = U0( I ) + WR( IR ) * ( U0( I+1 ) - U0( I ) )
               UT2 = U1( I ) + WR( IR ) * ( U1( I+1 ) - U1( I ) )

C              ------ Linear interpolation in time
               U = UT1 + WT * ( UT2 - UT1 )

               RTSRD( IR, ITOUT ) = RTSRD( IR, ITOUT ) + CONST * U
            END DO

         ELSE IF ( TOPOPT(4:4) .EQ. 'R' .AND. TIME .GE. TOUT(1) ) THEN

C           *** Case of an RTS (horizontal array) ***

C           ------ Linear interpolation in depth
            I = IRD( 1 )
            T = ( RD( 1 ) - Z( I ) ) / ( Z( I + 1 ) - Z( I ) )

            UT1 = U0( I ) + T * ( U0( I+1 ) - U0( I ) )
            UT2 = U1( I ) + T * ( U1( I+1 ) - U1( I ) )

C           ------ Linear interpolation in time
            U = UT1 + WT * ( UT2 - UT1 )

            CONST = SQRT( 2.0 ) * EXP( CI * RKT * V * TIME )

            DO IRNG = 1, NRR
               RM = 1000.0 * RR( IRNG )     ! Range in meters
C              REST = -SQRT(2.0) * COS( RKT*RM-PI/4.0 ) * SQRT(RKT/RM)
               REST = CONST * EXP( CI * ( -RKT * RM + PI / 4.0 ) ) *
     &                              SQRT( RKT / RM )
               RTSRR( IRNG, ITOUT ) =
     &            RTSRR( IRNG, ITOUT ) + DELTAK * U * REST
            END DO

         ENDIF
         ITOUT = ITOUT + 1

      END DO

      RETURN
      END
