C **********************************************************************
C *** SUBROUTINE REFLECT ***********************************************
C **********************************************************************
      SUBROUTINE REFL1(C3,C3S,B3,B3S)

      REAL IMAG

      DOUBLE PRECISION COSINE, ARGD, TWOPI

      COMPLEX E, G1, G2, SGAM2, ST2, V, Z1, ZT, CC3, CC3S, 
     &        CGAM1, CGAM2, SGAM3, CT1, Z, PHISIN


C*THIS SUBROUTINE CALCULATES REFLECTIVITY OF THE OCEAN BOTTOM AS A
C*FUNCTION OF THE GRAZING ANGLE.  THE BOTTOM IS A SMOOTH SOLID.
C  THE INPUT VALUES ARE:
C  C    : TOP-BOTTOM VELOCITY AT THE INTERFACE WITH THE BOTTOM
C  RHO3 : BOTTOM DENSITY
C  C3   : COMPRESSIONAL VELOCITY (BOTTOM SOUND VELOCITY)
C  CC3  : COMPLEX COMPRESSIONAL VELOCITY.
C  C3S  : SHEAR VELOCITY
C  CC3S : COMPLEX SHEAR VELOCITY.
C
      TWOPI=2.0D0*DACOS(-1.0D0)
C
C  CALCULATION OF COMPLEX COMPRESSIONAL VELOCITY.
C
      ALPHA=B3/(TWOPI*8.68588)
      IMAG=-ALPHA*C3
      CC3=CMPLX(C3,IMAG)
C
C  CALCULATION OF COMPLEX SHEAR SPEED
C
      ALPHA=B3S/(TWOPI*8.68588)
      IMAG=-ALPHA*C3S
      CC3S=CMPLX(C3S,IMAG)
      RETURN
C
C   THE REFLECTIVITY IS CALCULATED HERE.
C
      ENTRY REFL2(C,RHO2,RHO3,COSINE,RTHETA)
      ARGD=1.0D0-COSINE*COSINE
      ARGS=SIGN(1.0D0,ARGD)*MAX(ABS(ARGD),1.0D-30)
      PHISIN=CSQRT(CMPLX(ARGS,0.0))
      CGAM1=((CC3S/C)*SNGL(COSINE))**2
      G1=1.0-CGAM1
      SGAM2=CSQRT(G1)
      CT1=((CC3/C)*SNGL(COSINE))**2
      G2=1.0-CT1
      ST2=CSQRT(G2)
      Z=(RHO2*C)/PHISIN
      Z1=(RHO3*CC3)/ST2
      ZT=(RHO3*CC3S)/SGAM2
      CGAM2=(1.0-2.0*CGAM1)**2
      SGAM3=1.0-CGAM2
      E=Z1*CGAM2+ZT*SGAM3
      V=(E-Z)/(E+Z)
      IF(COSINE.LE.1.0)   RTHETA=(1.0-CABS(V)*CABS(V))
      IF(COSINE.GT.1.0)   RTHETA=2.0*AIMAG(V)
      RETURN
      END