      SUBROUTINE BCIMP( X, BCTYPE, BOTTOP, CPHS, CSHS, RHOHS,
     &   F, G, IPOW )

C     Compute Boundary Condition IMPedance

      INCLUDE 'comc.f'

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

      IPOW = 0

C     *** Get RHO, C just INSide the boundary ***
C     (There is at least one acoustic layer in the problem, except
C     in the case where BOUNCE is used to get the refl. coef. for
C     a purely elastic stack of layers.

      IF ( BOTTOP(1:3) .EQ. 'TOP' ) THEN
        IF ( NFACT .GT. 0 ) THEN
          I = LOC( NFACT ) + N( NFACT ) + 1
          RHOINS = RHO( I )
          CINS = PEKRT( OMEGA2 * H( NFACT ) **2 /
     &                  ( 2.0 + B1( NFACT ) ) )
        ENDIF
      ELSE
        IF ( NLACT .GT. 0 ) THEN
          I = LOC( NLACT ) + N( NLACT ) + 1
          RHOINS = RHO( I )
          CINS = PEKRT( OMEGA2 * H( NLACT ) **2 /
     &                  ( 2.0 + B1( NLACT ) ) )
        ENDIF
      ENDIF

C     *** Vacuum ***

      IF ( BCTYPE(1:1) .EQ. 'V' ) THEN
         F = 1.0
C         G = -CI * PEKRT( OMEGA2 / CINS ** 2 - X ) * SIGMA( 1 ) ** 2
         G = 0.0
      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 = PEKRT( GAMS2 )
            GAMP = PEKRT( 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 = PEKRT( 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        --- Shoot 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        --- Shoot 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 'comc.f'

      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 'comc.f'

      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
