      PROGRAM SCOOTER

C     Finite element fast-field program
C     Michael B. Porter   8/85

      INCLUDE 'COMSCO.FOR'

      COMPLEX GREEN( 300000 )

      CALL CLTIME
      CALL GETPAR

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

      CALL INIT
      CALL KERNEL( GREEN )
      CLOSE( GRNFIL )

      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 'COMSCO.FOR'

      CHARACTER*80 TITLE
      REAL*4 ZMIN, ZMAX

      TITLE = 'SCOOTER- '
      CALL READIN( TITLE, FREQ, MAXMED, NMEDIA,
     &   TOPOPT, CPT, CST, RHOT, BUMDEN, ETA, XI, N, SIGMA, DEPTH, 
     &   BOTOPT, CPB, CSB, RHOB,
     &   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

      IF ( RMAX .LE. 0.0 ) THEN
         CALL ERROUT( PRTFIL, 'F', 'GETPAR',
     &      'RMAX must be positive'  )
      ENDIF

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

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

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

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

C     *** Integration parameters ***

      RKMIN = 0.0
      IF ( CHIGH .NE. 0.0   ) RKMIN = SQRT( OMEGA2 ) / CHIGH
      IF ( CHIGH .GE. 1.0E6 ) RKMIN = 0.0

      IF ( CLOW .LE. 0.0 ) THEN
         WRITE( PRTFIL, * ) 'CLOW = ', CLOW
         CALL ERROUT( PRTFIL, 'F', 'GETPAR',
     &      'CLOW must be positive'  )
      ENDIF

      RKMAX = SQRT( OMEGA2 ) / CLOW

      NKPTS = 1000.0 * RMAX * ( RKMAX - RKMIN ) / PI
      WRITE( PRTFIL, * ) 'NKPTS = ', NKPTS

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

C     --- Set-up the vector of k-space points
      DELTAK = ( RKMAX - RKMIN ) / ( NKPTS - 1 )
      ATTEN  = DELTAK
      DO IK = 1, NKPTS
         RK( IK ) = RKMIN + ( IK - 1 ) * DELTAK
      END DO

      NSR = NSD * NRD
      NKCORE = MIN( 300000 / NSR, NKPTS, 4094 )
C     --- Write header for Green's function file
      CALL HEADER( TITLE, DELTAK, FREQ )

C     Optionally read in bottom reflection coefficient data

      DEGRAD = PI / 180.0
      IF ( BOTOPT(1:1) .EQ. 'F' ) THEN
         WRITE( PRTFIL, * )
     &                   'Using tabulated bottom reflection coefficient'
         READ( BRCFIL, * )
     &      NBPTS, ( THETAB( I ), RB( I ), PHIB( I ), I = 1, NBPTS )
         DO I = 1, NBPTS
            PHIB(I) = DEGRAD * PHIB(I)
         END DO
      ENDIF

C     Optionally read in top reflection coefficient data

      IF ( TOPOPT(2:2) .EQ. 'F' ) THEN
         WRITE( PRTFIL, * ) 'Using tabulated top reflection coefficient'
         READ( TRCFIL, * )
     &      NTPTS, ( THETAT( I ), RT( I ), PHIT( I ), I = 1, NTPTS )
         DO I = 1, NTPTS
            PHIT( I ) = DEGRAD * PHIT( I )
         END DO
      ENDIF

C     Optionally read in internal reflection coefficient data

      IF ( BOTOPT(1:1) .EQ. 'P' ) THEN
         WRITE( PRTFIL, * ) 'Reading precalculated refl. coeff. table'
         READ( IRCFIL ) NREF,
     &      ( XTAB(IK), FTAB(IK), GTAB(IK), ITAB(IK), IK = 1, NREF )
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE HEADER( TITLE, DELTAK, FREQ )

C     Routine to write header to disk file

      INCLUDE 'COMSCO.FOR'

      CHARACTER TITLE*80, PLTTYP*10

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

      OPEN ( FILE = 'GRNFIL', UNIT = GRNFIL, STATUS = 'UNKNOWN',
     &   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 ), NSD, NRD, NKPTS,
     &   NKCORE, REAL( DELTAK ), REAL( ATTEN ), NRECSR
      WRITE( GRNFIL, REC = 4 ) ( SD( I ), I = 1, NSD )
      WRITE( GRNFIL, REC = 5 ) ( RD( I ), I = 1, NRD )

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

      RETURN
      END
C**********************************************************************C
      SUBROUTINE INIT

C     Initializes arrays defining difference equations

      INCLUDE 'COMSCO.FOR'

      COMPLEX*16 CP( MAXN ), CS( MAXN ), CP2, CS2

      COMMON /SPACE/ CP, CS

      CMIN     = 1.0E6
      NFACT    = 0
      LOC( 1 ) = 0

C     *** Loop over media ***

      DO MED = 1, NMEDIA

         IF ( MED .GT. 1 .AND. SIGMA( MED ) .NE. 0.0 ) THEN
            WRITE( PRTFIL, * )
     &           'FATAL ERROR: Rough interfaces not allowed'
            STOP 'FATAL ERROR: Rough interfaces not allowed'
         ENDIF

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

         N1 = N(   MED ) + 1
         I  = LOC( MED ) + 1

         IF ( LOC( MED ) + N1 .GT. MAXN ) THEN
            WRITE( PRTFIL, * )
     &           'FATAL ERROR: Insufficient storage for mesh'
            STOP 'FATAL ERROR: Insufficient storage for mesh'
         ENDIF

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

         IF ( CS( I ) .EQ. ( 0.0, 0.0 ) ) THEN

C           *** Case of an acoustic medium ***

            MATER( MED )  =  'ACOUSTIC'
            IF ( NFACT .EQ. 0 ) NFACT = MED
            NLACT = MED

            DO J = I, I + N( MED )
               CMIN = MIN( CMIN, DBLE( CP( J ) ) )
               B1( J ) = OMEGA2 / CP( J ) ** 2
            END DO
         ELSE

C           *** Case of an elastic medium ***

            MATER( MED ) = 'ELASTIC'
            TWOH = 2.0 * H( MED )

            DO J = I, I + N( MED )
               CMIN = MIN( DBLE( CS( J ) ), CMIN )

               CP2 = CP( J ) ** 2
               CS2 = CS( J ) ** 2

               B1( J ) = TWOH / ( RHO( J ) * CS2 )
               B2( J ) = TWOH / ( RHO( J ) * CP2 )
               B3( J ) = 4.0 * TWOH * RHO(J) * CS2 * ( CP2 - CS2 ) / CP2
               B4( J ) = TWOH * ( CP2 - 2.0 * CS2 ) / CP2

               RHO( J ) = TWOH * OMEGA2 * RHO( J )
            END DO

         ENDIF
      END DO   ! next MED

      RETURN
      END
C**********************************************************************C
      SUBROUTINE BCIMP( X, BCTYPE, BOTTOP, CPHS, CSHS, RHOHS,
     &   F, G, IPOW )

C     Compute Boundary Condition IMPedance
C     Same subroutine as in KRAKENC except
C        PEKRT is replaced by SQRT
C        COMC  is replaced by COMSCO
C        CINS  is related to B1 differently

      INCLUDE 'COMSCO.FOR'

      COMPLEX*16 X, F, G, KX, KZ, TWERSK,
     &   GAMS2, GAMP2, GAMS, GAMP, RMU, YV( 5 ), RCMPLX
      COMPLEX*16 CPHS, CSHS, CINS
      CHARACTER BCTYPE*1, BOTTOP*3

      IPOW = 0

C     *** Get RHO, C just INSide the boundary ***

      IF ( BOTTOP(1:3) .EQ. 'TOP' ) THEN
         I = 1
         RHOINS = RHO( I )
         CINS   = SQRT( OMEGA2 / B1( 1 ) )
      ELSE
         I = LOC( NLACT ) + N( NLACT ) + 1
         RHOINS = RHO( I )
         CINS   = SQRT( OMEGA2 / B1( NLACT ) )
      ENDIF

C     *** Vacuum with Kirchoff roughness ***

      IF ( BCTYPE(1:1) .EQ. 'V' ) THEN
         F = 1.0
         G = -CI * SQRT( OMEGA2 / CINS ** 2 - X ) * SIGMA( 1 ) ** 2
      ENDIF

C     *** Vacuum with Twersky scatter model ***

      IF ( BCTYPE(1:1) .EQ. 'S' .OR. BCTYPE(1:1) .EQ. 'H' .OR.
     &     BCTYPE(1:1) .EQ. 'T' .OR. BCTYPE(1:1) .EQ. 'I' ) THEN

         OMEGA = SQRT( OMEGA2 )
         KX = SQRT( X )
         F = 1.0
         G = TWERSK( BCTYPE, OMEGA, BUMDEN, XI, ETA, KX, RHOINS, CINS )
         G = G / ( CI * OMEGA * RHOINS )
      ENDIF

C     *** Rigid ***
      IF ( BCTYPE(1:1) .EQ. 'R' ) THEN
         F = 0.0
         G = 1.0
      ENDIF

      YV(1) = F
      YV(2) = G
      YV(3) = 0.0
      YV(4) = 0.0
      YV(5) = 0.0

C     *** Acousto-elastic half-space ***

      IF ( BCTYPE(1:1) .EQ. 'A' ) THEN
         IF ( REAL( CSHS ) .GT. 0.0 ) THEN
            GAMS2 = X - OMEGA2 / CSHS ** 2
            GAMP2 = X - OMEGA2 / CPHS ** 2
            GAMS = SQRT( GAMS2 )
            GAMP = SQRT( GAMP2 )
            RMU = RHOHS * CSHS ** 2

            YV(1) = ( GAMS*GAMP - X ) / RMU
            YV(2) = ( ( GAMS2 + X )**2 - 4.0*GAMS*GAMP*X )*RMU
            YV(3) = 2.0*GAMS*GAMP - GAMS2 - X
            YV(4) = GAMP * ( X - GAMS2 )
            YV(5) = GAMS * ( GAMS2 - X )

            F = OMEGA2*YV( 4 )
            G = YV( 2 )
         ELSE
            GAMP = SQRT( X - OMEGA2 / CPHS**2 )
            F = 1.0
            G = RHOHS / GAMP
         ENDIF
      ENDIF

C     *** Tabulated reflection coefficient ***

      IF ( BCTYPE(1:1) .EQ. 'F' ) THEN
C        --- Compute the grazing angle THETA
         KX = SQRT( X )
         KZ = SQRT( OMEGA2 / CINS**2 - KX**2 )
         RADDEG = 360.0 / ( 2.0 * PI )
         THEINT = RADDEG * DATAN2( DBLE( KZ ), DBLE( KX ) )

C        --- Evaluate R(THEINT)
         IF ( BOTTOP(1:3) .EQ. 'TOP' ) THEN
            CALL REFCO( THEINT, RINT, PHIINT, THETAT, RT, PHIT, NTPTS )
         ELSE
            CALL REFCO( THEINT, RINT, PHIINT, THETAB, RB, PHIB, NBPTS )
         ENDIF

C        --- Convert R(THETA) to (f,g) in Robin BC
         RCMPLX = RINT * EXP( CI * PHIINT )
         F = 1.0
         G = ( 1.0 + RCMPLX ) / ( CI * KZ * ( 1.0 - RCMPLX ) )
      ENDIF

C     *** A top BC has the sign flipped relative to a bottom BC ***

      IF ( BOTTOP(1:3) .EQ. 'TOP' ) G = -G

C     *** Precalculated reflection coef ***

      IF ( BCTYPE(1:1) .EQ. 'P' ) THEN
         CALL IRCINT( X, F, G, IPOW, XTAB, FTAB, GTAB, ITAB, NKPTS )
      ENDIF

C     *** Shoot through elastic layers ***

      IF ( BOTTOP(1:3) .EQ. 'TOP' ) THEN

C        --- Shooting down from top
         IF ( NFACT. GT. 1 ) THEN

            DO MED = 1, NFACT-1
               CALL ELASDN( X, YV, IPOW, MED )
            END DO

            F = OMEGA2 * YV( 4 )
            G = YV( 2 )
         ENDIF
      ELSE

C        --- Shooting up from bottom
         IF ( NLACT. LT. NMEDIA ) THEN

            DO MED = NMEDIA, NLACT+1, -1
               CALL ELASUP( X, YV, IPOW, MED )
            END DO

            F = OMEGA2 * YV( 4 )
            G = YV( 2 )
         ENDIF
      ENDIF

      RETURN
      END
C**********************************************************************C
      SUBROUTINE ELASUP( X, YV, IPOW, MED )

C     Propagates through an elastic layer using 
C     compound matrix formulation

      INCLUDE 'COMSCO.FOR'

      COMPLEX*16 X, XV(5), YV(5), ZV(5), TWOX, XB3, FOURHX
      PARAMETER ( ROOF = 1.0E5, FLOOR = 1.0E-5, IPOWR = 5, IPOWF = -5 )

C     Euler's method for first step

      TWOX = 2.0*X
      TWOH = 2.0*H(MED)
      FOURHX = 4.0*H(MED)*X
      J = LOC(MED) + N(MED) + 1

      XB3 = X*B3(J) - RHO(J)

      ZV(1) = YV(1) - 0.5*(   B1( J ) * YV( 4 ) - B2( J ) * YV( 5 ) )
      ZV(2) = YV(2) - 0.5*( -RHO( J ) * YV( 4 ) -     XB3 * YV( 5 ) )
      ZV(3) = YV(3) - 0.5*(      TWOH * YV( 4 ) + B4( J ) * YV( 5 ) )
      ZV(4) = YV(4) - 0.5*(   XB3*YV(1) + B2(J)*YV(2) -TWOX*B4(J)*YV(3))
      ZV(5) = YV(5) - 0.5*(RHO(J)*YV(1) - B1(J)*YV(2) -    FOURHX*YV(3))

C     Modified midpoint method

      DO I = N( MED ), 1, -1
         J = J-1

         DO L = 1, 5
            XV(L) = YV(L)
            YV(L) = ZV(L)
         END DO

         XB3 = X * B3( J ) - RHO( J )

         ZV(1) = XV(1) - (   B1( J ) * YV( 4 ) - B2( J) * YV( 5 ) )
         ZV(2) = XV(2) - ( -RHO( J ) * YV( 4 ) -    XB3 * YV( 5 ) )
         ZV(3) = XV(3) - ( TWOH * YV( 4 ) + B4( J ) * YV( 5 ) )
         ZV(4) = XV(4) - (  XB3 *YV(1) + B2(J)*YV(2) - TWOX*B4(J)*YV(3))
         ZV(5) = XV(5) - (RHO(J)*YV(1) - B1(J)*YV(2) -     FOURHX*YV(3))

C        Scale if necessary

         IF ( I .NE. 1 ) THEN

            IF ( ABS( DBLE( ZV( 2 ) ) ) .LT. FLOOR ) THEN
               DO L = 1, 5
                  ZV(L) = ROOF * ZV( L )
                  YV(L) = ROOF * YV( L )
               END DO
               IPOW = IPOW-IPOWR
            ENDIF

            IF ( ABS( DBLE( ZV( 2 ) ) ) .GE. ROOF ) THEN
               DO L = 1, 5
                  ZV(L) = FLOOR * ZV( L )
                  YV(L) = FLOOR * YV( L )
               END DO
               IPOW = IPOW-IPOWF
            ENDIF
         ENDIF
      END DO

C     Apply the standard filter at the terminal point

      DO J = 1, 5
         YV( J ) = ( XV( J ) + 2.0*YV( J ) + ZV( J ) ) / 4.0
      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE ELASDN( X, YV, IPOW, MED )

C     Propagates through an elastic layer using 
C     compound matrix formulation

      INCLUDE 'COMSCO.FOR'

      COMPLEX*16 X, XV(5), YV(5), ZV(5), TWOX, XB3, FOURHX
      PARAMETER ( ROOF = 1.0E5, FLOOR = 1.0E-5, IPOWR = 5, IPOWF = -5 )

C     Euler's method for first step

      TWOX   = 2.0 * X
      TWOH   = 2.0 * H( MED )
      FOURHX = 4.0 * H( MED ) * X
      J = LOC(MED)+1

      XB3 = X * B3(J) - RHO(J)

      ZV(1) = YV(1) + 0.5*(   B1( J ) * YV( 4 ) - B2( J ) * YV( 5 ) )
      ZV(2) = YV(2) + 0.5*( -RHO( J ) * YV( 4 ) -     XB3 * YV( 5 ) )
      ZV(3) = YV(3) + 0.5*(      TWOH * YV( 4 ) + B4( J ) * YV( 5 ) )
      ZV(4) = YV(4) + 0.5*(   XB3*YV(1) + B2(J)*YV(2) -TWOX*B4(J)*YV(3))
      ZV(5) = YV(5) + 0.5*(RHO(J)*YV(1) - B1(J)*YV(2) - FOURHX*YV(3))

C     Modified midpoint method

      DO I = 1, N( MED )
         J = J+1

         DO L = 1, 5
            XV(L) = YV(L)
            YV(L) = ZV(L)
         END DO

         XB3 = X*B3(J) - RHO(J)

         ZV(1) = XV(1) + (   B1( J ) * YV( 4 ) - B2( J ) * YV( 5 ) )
         ZV(2) = XV(2) + ( -RHO( J ) * YV( 4 ) -     XB3 * YV( 5 ) )
         ZV(3) = XV(3) + (      TWOH * YV( 4 ) + B4( J ) * YV( 5 ) )
         ZV(4) = XV(4) + (   XB3*YV(1) + B2(J)*YV(2) -TWOX*B4(J)*YV(3))
         ZV(5) = XV(5) + (RHO(J)*YV(1) - B1(J)*YV(2) -    FOURHX*YV(3))

C        Scale if necessary

         IF ( I .NE. N(MED) ) THEN

            IF ( ABS( DBLE( ZV( 2 ) ) ) .LT. FLOOR ) THEN
               DO L = 1, 5
                  ZV(L) = ROOF * ZV( L )
                  YV(L) = ROOF * YV( L )
               END DO
               IPOW = IPOW-IPOWR
            ENDIF

            IF ( ABS( DBLE( ZV( 2 ) ) ) .GE. ROOF ) THEN
               DO L = 1, 5
                  ZV(L) = FLOOR * ZV( L )
                  YV(L) = FLOOR * YV( L )
               END DO
               IPOW = IPOW - IPOWF
            ENDIF

         ENDIF
      END DO

C     Apply the standard filter at the terminal point

      DO J = 1, 5
         YV( J ) = ( XV( J ) + 2.0 * YV( J ) + ZV( J ) ) / 4.0
      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE KERNEL( GREEN )

C     Solve system for a sequence of k-values

      INCLUDE 'COMSCO.FOR'

      INTEGER ISD( MAXNSD ), IRD( MAXNRD )
      REAL Z( MAXN ), RHOELT( MAXN ), WS( MAXNSD ), WR( MAXNRD )
      COMPLEX GREEN( NSD, NRD, * ), DF( MAXN ), EF( MAXN ), BELT
      COMPLEX*16 X

      IREC = 5
      NTOT = 0

C     --- Tabulate z-coordinates
      J = 1
      Z( 1 ) = DEPTH( NFACT )

      DO MED = NFACT, NLACT

         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 )

C     *** Main loop: solve system repeatedly ***

      IKCORE = 0
      IFLUSH = 0   ! number of records already flushed

C     --- NRECSR is number of records per S/R combo
      NRECSR = ( NKPTS - 1 ) / NKCORE + 1

C     *** Assemble matrix ***

      J = 1
      L = LOC( NFACT ) + 1
      DF( 1 ) = 0.0

      DO MED = NFACT, NLACT
         DO I = 1, N( MED )
            RHOELT( L ) = ( RHO( L ) + RHO( L + 1 ) ) / 2.0
            RHOH = RHOELT( L ) * H( MED )
            BELT = H( MED ) * ( ( B1( L ) + B1( L + 1 ) ) / 2.0 )
     &                          / ( 12.0 * RHOELT( L ) )

            DF( J     ) = DF( J ) - 1.0 / RHOH + 5.0 * BELT
            DF( J + 1 ) =         - 1.0 / RHOH + 5.0 * BELT
            EF( J + 1 ) =           1.0 / RHOH +       BELT

            J = J + 1
            L = L + 1
         END DO
         L = L + 1
      END DO

C     *** Step through each point in k-space ***

      DO IK = 1, NKPTS
C        WRITE( 6, * ) 'IK, NKPTS', IK, NKPTS
         IKCORE = IKCORE + 1
         X = ( RK( IK ) + CI * ATTEN ) ** 2

C        *** Solve for G(k) ***

         CALL SOLVE( NTOT1, IS, WS, ISD, WR, IRD, X, GREEN, IKCORE,
     &      DF, EF, RHOELT )

C        --- Flush buffer if necessary
         IF ( IKCORE .EQ. NKCORE .OR. IK .EQ. NKPTS ) THEN
            DO IS = 1, NSD

               DO IR = 1, NRD
C                 --- Index of source/receiver combo
                  ISR = ( IS - 1 ) * NRD + IR
                  IREC = 6 + ISR * NRECSR + IFLUSH
                  WRITE( GRNFIL, REC = IREC )
     &               ( GREEN( IS, IR, IKCORE ), IKCORE = 1, NKCORE )
               END DO
 
            END DO
            IKCORE = 0
            IFLUSH = IFLUSH + 1
         ENDIF

      END DO

      RETURN
      END
C**********************************************************************C
      SUBROUTINE SOLVE( NTOT1, IS, WS, ISD, WR, IRD, X,
     &   GREEN, IKCORE, DF, EF, RHOELT )

C     Set up the linear system and solve

      INCLUDE 'COMSCO.FOR'

      INTEGER ISD( * ), IRD( * )
      REAL RHOELT( * ), WS( * ), WR( * )
      COMPLEX GREEN( NSD, NRD, * ), D( MAXN ), E( MAXN ),
     &   B( MAXN ), RV1( MAXN ), RV2( MAXN ), RV3( MAXN ), RV4( MAXN ),
     &   DF( * ), EF( * ), BELT, XT
      COMPLEX*16 X, F, G
      CHARACTER BCTYPE*1

      COMMON /SPACE/ RV1, RV2, RV3, RV4

      EQUIVALENCE ( D, B )

C     *** Assemble matrix ***

      J = 1
      L = LOC( NFACT ) + 1
      D( 1 ) = DF( 1 )

      DO MED = NFACT, NLACT
         XT = -H( MED ) * X / 12.0

         DO I = 1, N( MED )
            BELT = XT / RHOELT( L )
            D( J    ) = D(  J     ) + 5.0 * BELT
            D( J + 1) = DF( J + 1 ) + 5.0 * BELT
            E( J + 1) = EF( J + 1 ) +       BELT
            J = J + 1
            L = L + 1
         END DO

         L = L + 1
      END DO

C     *** Corner elt requires top impedance ***

      BCTYPE(1:1) = TOPOPT(2:2)
      CALL BCIMP( X, BCTYPE, 'TOP', CPT, CST, RHOT, F, G, IPOW )

      IF ( G .EQ. 0.0 ) THEN
         D( 1 ) = 1.0
         E( 2 ) = 0.0
      ELSE
         D( 1 ) = D( 1 ) + F / G
      ENDIF

C     *** Corner elt requires bottom impedance ***

      BCTYPE(1:1) = BOTOPT(1:1)
      CALL BCIMP( X, BCTYPE, 'BOT', CPB, CSB, RHOB, F, G, IPOW )

      IF ( G .EQ. 0.0 ) THEN
         D( NTOT1 ) = 1.0
         E( NTOT1 ) = 0.0
      ELSE
         D( NTOT1 ) =  D( NTOT1 ) - F / G
      ENDIF

C     * Do LU decomposition *

      CALL FACTOR( NTOT1, D, E, IERR, RV1, RV2, RV3, RV4 )

C     *** Loop over all source positions ***

      DO IS = 1, NSD

C        --- Set up RHS
         DO J = 1, NTOT1
            B( J ) = 0.0
         END DO

         I = ISD( IS )

C        --- depending on the normalization, we may need rho(zs)
         RHOSD = 1.0

         B( I     ) = 2.0 * ( 1.0 - WS( IS ) ) / RHOSD
         B( I + 1 ) = 2.0 *     WS( IS )       / RHOSD

C        * Solve the system and extract soln for desired receivers *

         CALL BACKSB( NTOT1, D, E, IERR, RV1, RV2, RV3, RV4, B )

         DO IR = 1, NRD
            I = IRD( IR )
            GREEN( IS, IR, IKCORE ) =
     &                B( I ) + WR( IR ) * ( B( I + 1 ) - B( I ) )
         END DO

      END DO

      RETURN
      END
