C 
C *** AMPLITUDE FACTER OF MULTIPLE REFLECTION METHOD OF SH WAVE
C
      PROGRAM MSHML2
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 IN SUBROUTINE !!!
C ** DNS(I) : DENSITY OF Ith LAYER
C ** VS (I) : VS OF Ith LAYER
C ** DEP(I) : DEPTH ( THICKNESS ) OF Ith LAYER
C ** NL     : TOTAL NUMBER OF LAYERS
C ** ML     : NL IN MAIN PROGRAM DIMENSION
C ** THETA  : INCIDENT ANGLE OF SH WAVE ( RADIAN )
C OUTPUT DATA  
C ** AMP(I) : AMPLITUDE FACTER ON SURFACE OF Ith OMEGA ( (I-1)*DOM+FOM )
C ** SPCA(I): SPECTRUM AMPLITUDE OF AMP(I)
C
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(ML=100,MOM=5000)
      DIMENSION DNS(ML),DEP(ML),VS(ML),QF(ML)
      COMPLEX*16 AMP(MOM)
      character name*80,fname*30
C
      WRITE(*,10)
   10 FORMAT(' Enter Soil Data File Name ')
      READ(*,'(A)')FNAME
      WRITE(*,20)
   20 FORMAT(' Enter Soil Data File Number ')
      READ(*,*)NNN      

      OPEN (4,FILE=FNAME)
      OPEN(7,FILE='spec_dat')

      do 100 iiii = 1, nnn

      CALL SOILIN(FOM,DOM,NOM,DNS,VS,DEP,QF,NQF,NL,ML,THETA,NSF,iiii,
     *            NAME)
      CALL SBSHML(FOM,DOM,NOM,MOM,DNS,VS,DEP,QF,NQF,NL,ML,THETA,NSF,AMP)
      CALL SPCOUT(FOM,DOM,NOM,MOM,AMP,NAME)

  100 continue

      CLOSE(4)
      CLOSE(7)
C
      STOP
      END
c
C
C *** READ OF SOIL DATA FROM DATA FILE
C 
      SUBROUTINE SOILIN(FOM,DOM,NOM,DNS,VS,DEP,QF,NQF,NL,ML,THETA,
     *                  NSF,iiii,NAME)
C
C ** FOM    : FIRST CIRCULAR FREQUENCY ( OMEGA )
C ** DOM    : INCREMENT OF CIRCULAR FREQUENCY ( DELTA OMEGA )
C ** NOM    : TOTAL NUMBER OF CIRCULAR FREQUENCY
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 = 1/(2*H)
C ** NQF    : = 1 : QF ENTER,   = 0 : QF NOT NOT ENTER
C ** NL     : TOTAL NUMBER OF LAYERS
C ** ML     : NL IN MAIN PROGRAM DIMENSION
C ** THETA  : INCIDENT ANGLE OF SH WAVE ( DEGREE --> RADIAN )
C ** NSF    : =1:OUTPUT POINT IS ON SURFACE, =0:OUTPUT POINT IS IN SOIL
C
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER FNAME*30,name*80
      DIMENSION DNS(ML),DEP(ML),VS(ML),QF(ML)
C
C *** SOIL DATA INPUT ***
C 
      if(iiii.eq.1) then
        READ(4,1000)FOM,DOM,NOM
        READ(4,1010)THET1
        READ(4,1020)NSF
      end if
      READ(4,1025)NAME
      write(*,'(a)')name
      READ(4,1030)NL,NQF
      DO 100 IL=1,NL
      READ(4,1040)DNS(IL),VS(IL),DEP(IL),QF(IL)
  100 CONTINUE
 1000 FORMAT(2F10.0,I10)
 1010 FORMAT(F10.0)
 1020 FORMAT(I10)
 1025 FORMAT(a)
 1030 FORMAT(2I10)
 1040 FORMAT(4F10.0)
C
C *** CHANGE THETA FROM DEGREE TO RADIAN ***
C
      THETA=THET1*0.017453292520D0
C
C *** DATA ECHO ***
C
      WRITE(*,2000)
      WRITE(*,2010)FOM,DOM,NOM
      WRITE(*,2020)THET1,THETA
      WRITE(*,2030)NSF
      WRITE(*,2040)NL, NQF
      DO 200 IL=1,NL
      WRITE(*,2050)IL,DNS(IL),VS(IL),DEP(IL),QF(IL),0.5/QF(IL)
  200 CONTINUE
      WRITE(*,*)
 2000 FORMAT(/' FROM SUBROUTINE SOILIN : INPUT DATA ECHO')
 2010 FORMAT(' FOM,DOM,NOM    =',2F10.5,I5)
 2020 FORMAT(' THETA (DEG) =',F10.4,'  (RAD) =',F10.4)
 2030 FORMAT(' = 1:OUTPUT POINT IS ON SURFACE, = 0:IN SOIL : ',I3)
 2040 FORMAT(' TOTAL RAYLERS NUMBER =',I5,' QF IN ( =1:IN )=',I3)
 2050 FORMAT(' IL,DNS,VS,DEP,QF,H =',I3,F6.3,3F10.3,F10.6)
C       
      RETURN
      END
c
C
C *** CUL. OF AMPLITUDE FACTER BY MULTIPLE REFLECTION METHOD OF SH WAVE
C 
      SUBROUTINE SBSHML(FOM,DOM,NOM,MOM,DNS,VS,DEP,QF,NQF
     *                 ,NL,ML,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
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 ** NL     : TOTAL NUMBER OF LAYERS
C ** ML     : NL IN MAIN PROGRAM DIMENSION
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)
      DIMENSION DNS(ML),DEP(ML),VS(ML),QF(ML)
      COMPLEX*16 AMP(MOM),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 *** MOM CHECK ***
C
C     IF(MOM.NE.32800)THEN
C     WRITE(*,*)' FROM SUBROUTINE SBSHML :'
C     WRITE(*,*)' MOM.NE.32800 !!!!, THEN STOP !!!!'
C     STOP
C     END IF
C
C *** CAL. OF AMPLITUDE 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
      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)
          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)=SF/(R11-R21*DCMPLX(0.0D0,1.0D0)/WNSN)
C
  200 CONTINUE
C
      RETURN
      END
c
C
C *** OUTPUT OF SPECTRUM DATA ON SCREEN AND FILE DATA ***
C
      SUBROUTINE SPCOUT(FOM,DOM,NOM,MOM,SPC,NAME)
C
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 ** SPC(I) : SPECTRUM ( COMPLEX*8 )
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 SPC(MOM)
      character name*80
C
C *** MOM CHECK ***
C
C     IF(MOM.NE.32800)THEN
C     WRITE(*,*)' FROM SUBROUTINE SPCOUT :'
C     WRITE(*,*)' MOM.NE.32800 !!!!, THEN STOP !!!!'
C     STOP
C     END IF
C
C *** SPECTRUM OUTPUT ON SCREEN ***
C 
c      WRITE(*,1000)
c      WRITE(*,1010)
c      READ(*,*)NYES
c      IF(NYES.NE.0)GOTO 10
C
c      WRITE(*,1020)
c      IFOM=1
c      IF(FOM.EQ.0.0D0)THEN
c      IFOM=2
c      WRITE(*,1030)1,FOM,0.0D0,SPC(1),CDABS(SPC(1))
c      END IF
C
c      DO 100 IOM=IFOM,NOM
c      OM=DFLOAT(IOM-1)*DOM+FOM
c      PER=6.2831853070D0/OM
c      WRITE(*,1030)IOM,OM,PER,SPC(IOM),CDABS(SPC(IOM))
c  100 CONTINUE
C
c      WRITE(*,*)
c 1000 FORMAT(/' FROM SUBROUTINE SPCOUT : SPECTRUM OUTPUT')
c 1010 FORMAT(' SPECTRUM DATA OUTPUT ON SCREEN ? ( YES = 0 ) : ')
c 1020 FORMAT(/' NUM., OMEGA, PERIOD,  SPC(REAL),  SPC(IMAG),  SPC(ABS)')
c 1030 FORMAT(I4,2F8.3,3D12.4)  
C
C *** FILE OUTPUT TO FILE DATA ***
C
c   10 WRITE(*,2000)
c 2000 FORMAT(' FILE OUTPUT ? ( YES = 0 : spec_dat ) ')
c      READ(*,*)NYES
c      IF(NYES.NE.0)GOTO 20
      NYES=0 
C

      IFOM=2

      write(7,'(a)')name
      WRITE(7,2100)FOM,CDABS(SPC(1))
      DO 200 IOM=IFOM,NOM
      OM=DFLOAT(IOM-1)*DOM+FOM
      WRITE(7,2100)om,CDABS(SPC(IOM))
  200 CONTINUE
 2100 FORMAT(2F10.4)
   20 WRITE(*,*)
C
      RETURN
      END
