C
C *** CUL. OF AMPLITUDE FACTER BY MULTIPLE REFLECTION METHOD OF SH WAVE
C 
      SUBROUTINE SBMLTPL(FOM,DOM,NOM,MOM,THETA,NSF,AMP)
C
C INPUT DATA  
C ** FOM    : FIRST CIRCULAR FREQUENCY ( OMEGA )
C ** DOM    : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C ** NOM    : TOTAL NUMBER OF CIRCULAR FREQUENCY
C ** MOM    : NOM IN MAIN PROGRAM DIMENSION
C           : IF MOM > 2**13(=8192), CAN NOT USE !!!
C ** DNS(I) : DENSITY OF Ith LAYER              ==> COMMON
C ** VS (I) : VS OF Ith LAYER                   ==> COMMON
C ** DEP(I) : DEPTH ( THICKNESS ) OF Ith LAYER  ==> COMMON
C ** QF(I)  : QUALITY FACTER OF Ith LAYER       ==> NEGRECT
C ** NQF    : = 0 : QF IN,   = 1 : QF NOT IN    ==> NQF = 1
C ** NL     : TOTAL NUMBER OF LAYERS            ==> COMMON
C ** THETA  : INCIDENT ANGLE OF SH WAVE (RADIAN)
C ** NSF    : = 1:OUTPUT POINT IS ON SURFACE,  = 0:IN SOIL
C OUTPUT DATA  
C ** AMP(I) : AMPLITUDE FACTER ON SURFACE OF Ith OMEGA ( (I-1)*DOM+FOM )
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /LAYER/ NL,VS(35),DEP(35),SM(35) /DNSITY/ DNS(35)
      COMPLEX*16 AMP(MOM),VSN2,CINV2,R11,R12,R21,R22,VS2,SMC,WN,WND,WNS
     *      ,COWN,SIWN,SA11,SA12,SA21,SA22,SB11,SB12,SB21,SB22,WNN,WNSN
C
C *** CAL. OF AMPLITUDE FACTER ***
C
      VSN2=DCMPLX(VS(NL)*VS(NL),0.0D0)
C     IF(NQF.EQ.1)VSN2=VSN2*DCMPLX(1.0D0,1.0D0/QF(NL))
      TSIN=DSIN(THETA)
      CINV2=TSIN*TSIN/VSN2
C
      SF=1.0D0
      IF(NSF.EQ.1)SF=2.0D0
C
      IFOM=1
      IF(FOM.EQ.0.0D0)THEN
      AMP(1)=DCMPLX(SF,0.0D0)
      IFOM=2
      END IF
C
      DO 200 IO=IFOM,NOM
C
        OM=DFLOAT(IO-1)*DOM+FOM
C
        R11=(1.0D0,0.0D0)
        R12=(0.0D0,0.0D0)
        R21=(0.0D0,0.0D0)
        R22=(1.0D0,0.0D0)
C
        DO 100 IL=NL-1,1,-1
C
          SA11=R11
          SA12=R12
          SA21=R21
          SA22=R22
C
          VS2=DCMPLX(VS(IL)*VS(IL),0.0D0)
C         IF(NQF.EQ.1)VS2=VS2*DCMPLX(1.0D0,1.0D0/QF(IL))
          SMC=DNS(IL)*VS2   
          WN=OM*CDSQRT(1.0D0/VS2-CINV2)
          WND=WN*DEP(IL)
          WNS=WN*SMC
          COWN=CDCOS(WND)
          SIWN=CDSIN(WND)
C
          SB11=COWN
          SB12=SIWN/WNS
          SB21=-WNS*SIWN
          SB22=COWN
C
          R11=SA11*SB11+SA12*SB21
          R12=SA11*SB12+SA12*SB22
          R21=SA21*SB11+SA22*SB21
          R22=SA21*SB12+SA22*SB22
C
  100   CONTINUE
C
        WNN=OM*CDSQRT(1.0D0/VSN2-CINV2)
        WNSN=WNN*DNS(NL)*VSN2
C
        AMP(IO)=SF/(R11-R21*DCMPLX(0.0D0,1.0D0)/WNSN)
C
  200 CONTINUE
C
      RETURN
      END
