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

C     Compute Boundary Condition IMPedance

      INCLUDE 'comc1.f'

      COMPLEX*16 X, F, G, KX, KZ, TWERSK,
     &   GAMS2, GAMP2, GAMS, GAMP, RMU, YV( 5 ),
     &   PEKRTZ, RCMPLX
      COMPLEX*16 CPHS, CSHS, CINS
      REAL*8 RSX, RKS, RKP, OMEGA, CENT
      CHARACTER BCTYPE*1, BOTTOP*3

      IPOW = 0
      OMEGA = SQRT( OMEGA2 )

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

      IF ( BOTTOP(1:3) .EQ. 'TOP' ) THEN
         I = LOC( NFACT ) + N( NFACT ) + 1
         RHOINS = RHO( I )
         CINS = PEKRT( OMEGA2 * H( NFACT ) **2 / ( 2.0 + B1( NFACT )))
      ELSE
         I = LOC( NLACT ) + N( NLACT ) + 1
         RHOINS = RHO( I )
         CINS = PEKRT( OMEGA2 * H( NLACT ) **2 / ( 2.0 + B1( NLACT )))
      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
         RSX = DREAL( SQRT( X ) )
         IF ( REAL( CSHS ) .GT. 0.0 ) THEN
            RKS = OMEGA * DREAL( CSHS ) / CSHS / DCONJG( CSHS )
            RKP = OMEGA * DREAL( CPHS ) / CPHS / DCONJG( CPHS )
            GAMS2 = X - OMEGA2 / CSHS ** 2
            GAMP2 = X - OMEGA2 / CPHS ** 2

            GAMS = PEKRTZ( GAMS2, CENT, RKS )
           GAMP = PEKRTZ( GAMP2, CENT, RKP )
c            GAMS = PEKRT( GAMS2 )
c            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
            RKP = OMEGA * DREAL( CPHS ) / CPHS / DCONJG( CPHS )
c            write(*,*) ' RKP : ', RKP
            GAMP2 = X - OMEGA2 / CPHS ** 2
           GAMP = PEKRTZ( GAMP2, CENT, RKP )
c            GAMP = PEKRT( GAMP2 )        
              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
      