C
C *** CUL. OF AMPLITUDE FACTER BY MULTIPLE REFLECTION METHOD OF SH WAVE
C 
      SUBROUTINE WVML(FOM,DOM,NOM,NDLYR,NL,DNS,VS,DEP,QF,NQF
     *               ,THETA,NDFFT,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 NDLYR  : NL IN MAIN PROGRAM DIMENSION
C NL     : TOTAL NUMBER OF LAYERS
C DNS(I) : DENSITY OF Ith LAYER
C VS (I) : VS OF Ith LAYER
C DEP(I) : DEPTH (THICKNESS) OF Ith LAYER
C QF(I)  : QUALITY FACTER OF Ith LAYER
C NQF    : = 0 : QF IN,   = 1 : QF NOT IN
C THETA  : INCIDENT ANGLE OF SH WAVE (RADIAN)
C NDFFT  : 
C OUTPUT DATA  
C AMP(I) : AMPLIFICATION FACTER OF Ith OMEGA ((I-1)*DOM+FOM)
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DNS(NDLYR),DEP(NDLYR),VS(NDLYR),QF(NDLYR)
      COMPLEX*16 AMP(NDFFT),VSN2,CINV2,R11,R12,R21,R22,VS2,SM,WN,WND,WNS
     *      ,COWN,SIWN,SA11,SA12,SA21,SA22,SB11,SB12,SB21,SB22,WNN,WNSN
C
C *** CAL. OF AMPLIFICATION FACTER ***
C
      VSN2=DCMPLX(VS(NL)*VS(NL),0.0D0)
      IF(NQF.EQ.1)VSN2=VSN2*DCMPLX(1.0D0,1.0D0/QF(NL))
      TSIN=DSIN(THETA)
      CINV2=TSIN*TSIN/VSN2
C
      IFOM=1
      IF(FOM.EQ.0.0D0)THEN
      AMP(1)=DCMPLX(2.D0,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)
          IF(NQF.EQ.1)VS2=VS2*DCMPLX(1.0D0,1.0D0/QF(IL))
          SM=DNS(IL)*VS2   
          WN=OM*CDSQRT(1.0D0/VS2-CINV2)
          WND=WN*DEP(IL)
          WNS=WN*SM
          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)=2.D0/(R11-R21*DCMPLX(0.0D0,1.0D0)/WNSN)
C
  200 CONTINUE
C
      RETURN
      END
