      PROGRAM SHML
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION DM(100),AAMP(100)
      COMPLEX*16 AMP,SA11,SA12,SA21,SA22,SB11,SB12,SB21,SB22,
     *           R11,R12,R21,R22,G(100),VS2(100),KM(100),
     *           KD,COKD,SIKD,GKM
C
      PAI=DATAN(1.D0)*4.D0
C
C     CALL CLOCK
C
      WRITE(*,*) 'FILE-OUT  0 OR 1 (I1) ?'
      READ(*,5030) NFO
      WRITE(*,*) 'HOW MANY LAYERS (I2)  ?'
      READ(*,5000) NL
      WRITE(*,*) 'INCIDENT ANGLE  ?'
      READ(*,5010) THETA
      THETA=THETA*DATAN(1.D0)/45.D0
      IF(NFO.EQ.1) WRITE(7,7000) NL,THETA
      DO 100 J=1,NL
        WRITE(*,*) ' '
        WRITE(*,6000) J,'....DENSITY  ?'
        READ(*,5010) DENS
        WRITE(*,6000) J,'....VS       ?'
        READ(*,5010) VS
        WRITE(*,6000) J,'.... H       ?'
        READ(*,5010)  H
C       H=2.D0*PAI/QV
        G0=VS*VS*DENS
        G(J)=G0*(1.D0+2.D0*DCMPLX(0.D0,H))
        VS2(J)=G(J)/DENS
        IF(J.EQ.NL) THEN
          IF(NFO.EQ.1) WRITE(7,7010) DENS,VS,H
          GO TO 100
        END IF
        WRITE(*,6000) J,'....DEPTH    ?'
        READ(*,5010) DM(J)
        IF(NFO.EQ.1) WRITE(7,7010)DENS,VS,H,DM(J)
  100 CONTINUE
      WRITE(*,*) ' '
      WRITE(*,*) 'D OMEGA      ?'
      READ(*,5010) DO
      WRITE(*,*) 'NO  (I3)     ?'
      READ(*,5020) NO
      IF(NFO.EQ.1) WRITE(7,7020) DO,NO
C
      WRITE(*,6005)
C
      DO 900 I=1,NO
C
        OM=DFLOAT(I)*DO
C
C
        R11=DCMPLX(1.D0,0.D0)
        R12=DCMPLX(0.D0,0.D0)
        R21=DCMPLX(0.D0,0.D0)
        R22=DCMPLX(1.D0,0.D0)
C
        DO 300 J=NL-1,1,-1
C
          SA11=R11
          SA12=R12
          SA21=R21
          SA22=R22
          KM(J)=OM*CDSQRT(1.D0/VS2(J)-DSIN(THETA)*DSIN(THETA)/VS2(NL))
          KD=KM(J)*DM(J)
          COKD=CDCOS(KD)
          SIKD=CDSIN(KD)
          GKM=G(J)*KM(J)
          SB11=COKD
          SB12=SIKD/GKM
          SB21=-GKM*SIKD
          SB22=COKD
          R11=SA11*SB11+SA12*SB21
          R12=SA11*SB12+SA12*SB22
          R21=SA21*SB11+SA22*SB21
          R22=SA21*SB12+SA22*SB22
C
  300   CONTINUE
C
        KM(NL)=OM*CDSQRT((1.D0-DSIN(THETA)*DSIN(THETA))/VS2(NL))
C
        AMP=2.D0/(R11-R21*DCMPLX(0.D0,1.D0)/KM(NL)/G(NL))
C
        WRITE(*,6010) OM,2*PAI/OM,AMP,CDABS(AMP)
C
        IF(NFO.EQ.1) THEN
          AAMP(I)=CDABS(AMP)
        END IF
C
  900 CONTINUE
C
      IF(NFO.EQ.1) THEN
        DO 1000 I=1,NO,10
          WRITE(7,7030) AAMP(I),AAMP(I+1),AAMP(I+2),AAMP(I+3),
     *                  AAMP(I+4),AAMP(I+5),AAMP(I+6),AAMP(I+7),
     *                  AAMP(I+8),AAMP(I+9)
 1000   CONTINUE
      END IF
C
C
C    CALL CLOCK(TIME)
C
      WRITE(*,6020) TIME
C
 5000 FORMAT(I2)
 5010 FORMAT(F15.3)
 5020 FORMAT(I3)
 5030 FORMAT(I1)
 6000 FORMAT(I2,A15)
 6005 FORMAT(//' OM(RAD/SEC),PERIOD(SEC),AMP(REAL),AMP(IMAG),AMP(ABS)'/)
 6010 FORMAT(F7.3,F8.3,2F11.4,F12.4)
 6020 FORMAT(//5X,'CPU TIME ',F15.7,' (SEC)')
 7000 FORMAT(I15,F15.5)
 7010 FORMAT(2F8.2,F8.5,F8.2)
 7020 FORMAT(F15.7,I15)
 7030 FORMAT(10F8.2)
C
      STOP
      END
