      PROGRAM RAY2LE
C
C *** CAL. OF PHASE & GROUP VELOCITY, AMPLITUDE RATIO ( U0/W0 )
C     AND AMPLITUDE SPECTRUM OF RAYLEIGH WAVE IN TWO LAYERS MEDIUM ***
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CC(0:100),UU(0:100),CHECK(0:100)
     *       ,RKL(0:100),ARK(0:100)
C
C     CALL CLOCK
C
C *** INPUT DATA ***
C
    5 WRITE(*,*) 'FILE-OUT  0 OR 1 (I1) ?'
      READ(*,*)NFO
C       WRITE(*,*) ' '
        WRITE(*,*) '....DENSITY 1 ?'
        READ(*,*) DN1
        WRITE(*,*) '....VP 1     ?'
        READ(*,*) VP1
        WRITE(*,6000) '....VS 1     ?'
        READ(*,*) VS1
        WRITE(*,6000) '....DENSITY 2 ?'
        READ(*,*) DN2
        WRITE(*,6000) '....VP 2     ?'
        READ(*,*) VP2
        WRITE(*,6000) '....VS 2     ?'
        READ(*,*) VS2
        WRITE(*,6000) '....DEPTH    ?'
        READ(*,*) DM
C
      WRITE(*,*) 'FUNDAMENTAL MODE ? ( YES = 0 (II) !!! )'
      READ(*,*) NFUN
      IF(NFUN.NE.0)THEN
      WRITE(*,*) 'DO YOU HAVE INITIAL OM AND CC ? ( YES = 0 (II) !!! )'
      READ(*,*) NIC
      IF(NIC.NE.0)WRITE(*,*)'PLEASE USE RAY2LN !!! AND STOP !!!'
      IF(NIC.NE.0)STOP
      WRITE(*,*) 'INITIAL OM   ?'
      READ(*,*) FOM
      WRITE(*,*) 'INITIAL CC   ?'
      READ(*,*) FCC
      ELSE
      FOM=0.D0
      END IF
C
      WRITE(*,*) 'D OMEGA      ?'
      READ(*,*) DO
   10 WRITE(*,*) 'NO  (I3)     ?'
      READ(*,*) NO
      EOM=DO*NO+FOM
      IF(NFUN.NE.0)THEN
      WRITE(*,*) 'FINAL OMEGA =',EOM,' OK ? (YES= 0 (I1))'
      READ(*,*) NOK
      IF(NOK.NE.0)GOTO 10
      END IF
      WRITE(*,*) 'NUM.OF DDO (I3) ?'
      READ(*,*) NDO
      WRITE(*,*) 'DELTA C      ?'
      READ(*,*) DC
C
 6000 FORMAT(A15)
 5010 FORMAT(F10.0)
 5020 FORMAT(I3)
 5030 FORMAT(I1)
C
C     DN1=2.0D0
C     VP1=2000.0D0
C     VS1=1000.0D0
C     DN2=2.50D0
C     VP2=5000.0D0
C     VS2=3000.0D0
C     DM =2500.0D0
      DZ =0.010D0
C
      WRITE(*,1000)
      WRITE(*,1010)NFO
      WRITE(*,1015)NFUN
      WRITE(*,1020)VP1,VS1,DN1
      WRITE(*,1030)VP2,VS2,DN2
      WRITE(*,1040)DM,DZ,DC
      WRITE(*,1050)DO,NO,NDO
C
      P21=VP1*VP1
      S21=VS1*VS1
      P22=VP2*VP2
      S22=VS2*VS2
      PS1=(P21-2.D0*S21)/2.D0/(P21-S21)
      PS2=(P22-2.D0*S22)/2.D0/(P22-S22)
      WRITE(*,*) '1L POISSON''S RATIO =',PS1
      WRITE(*,*) '2L POISSON''S RATIO =',PS2
      WRITE(*,*)
      WRITE(*,*) 'ALL INPUT DATA OK ? (YES= 0 (I1))'
      READ(*,5020) NOK2
      IF(NOK2.NE.0)GOTO 5
C
      IF(NFO.EQ.1) THEN
        WRITE(7,7010) DN1,VP1,VS1,DM
        WRITE(7,7010) DN2,VP2,VS2
        WRITE(7,7020) DO,NO,NDO,NFUN
      END IF
C
C *** Cul. of CR2 ,CR1 and CC ***
C
      CR2=FNCCR(VP2,VS2,DZ)
      CR1=FNCCR(VP1,VS1,DZ)
C
      IF(NFUN.EQ.0)FCC=CR2
      C0=FCC
      CHECK(0)=FNCFR2(FOM,C0,VP1,VS1,DN1,VP2,VS2,DN2,DM)
      WRITE(*,1100)CR2,CR1,CHECK(0)
C
      CC(0)=FCC
      DO 100 I=1,NO*NDO+1
      OM=DFLOAT(I)*DO/DFLOAT(NDO)+FOM
      RF0=FNCFR2(OM,C0,VP1,VS1,DN1,VP2,VS2,DN2,DM)
  150 C1=C0-DC
      RF1=FNCFR2(OM,C1,VP1,VS1,DN1,VP2,VS2,DN2,DM)
      AF=RF0/RF1
      C0 =C1
      R00=RF0
      RF0=RF1
      IF(AF.GE.0.D0)GOTO 150
      CC(I)=C1+DC*DABS(RF0/(RF0-R00))
      C0=CC(I)
      CHECK(I)=FNCFR2(OM,C0,VP1,VS1,DN1,VP2,VS2,DN2,DM)
  100 CONTINUE
C
C *** Cal. of UU ***
C
      UU(0)=CR2
      OM=DO/10000.D0+FOM
      CALL UWU(OM,CC(0),VP1,VS1,DN1,VP2,VS2,DN2,DM,U1,R1,AK)
      RKL(0)=R1
      ARK(0)=AK
      DO 200 I=NDO,NO*NDO,NDO
      OM=DFLOAT(I)*DO/DFLOAT(NDO)+FOM
      CALL UWU(OM,CC(I),VP1,VS1,DN1,VP2,VS2,DN2,DM,U1,R1,AK)
      UU(I)=U1
      RKL(I)=R1
      ARK(I)=AK
  200 CONTINUE
C
C *** OUTPUT DATA ***
C
      WRITE(*,1200)
      DO 300 I=0,NO*NDO,2*NDO
      I1=I
      I2=I+NDO

      DFN=DO/DFLOAT(NDO)
      OM1=DFLOAT(I1)*DFN+FOM
      OM2=DFLOAT(I2)*DFN+FOM

      WRITE(*,1210)OM1,CC(I1),UU(I1),RKL(I1),ARK(I1)
     *            ,OM2,CC(I2),UU(I2),RKL(I2),ARK(I2)
  300 CONTINUE
C
      IF(NFO.EQ.1)THEN
      WRITE(7,7030)CR2,CR1
      IF(NFUN.NE.0)WRITE(7,7010)FOM,FCC
      DO 400 I=0,NO*NDO,2*NDO
      I1=I
      I2=I+NDO
      UK1=ARK(I1)*RKL(I1)
      UK2=ARK(I2)*RKL(I2)
      WRITE(7,7030)CC(I1),UU(I1),ARK(I1),UK1
     *            ,CC(I2),UU(I2),ARK(I2),UK2
  400 CONTINUE
      END IF
C
C     CALL CLOCK(TIME)
C     WRITE(*,2000)TIME
C
 1000 FORMAT(//' Rayreigh Wave C for 2 Layers'/)
 1010 FORMAT(' NFO( OUTPUT FILE ) =',I5)
 1015 FORMAT(' NFUN =',I2,': 0 (FUNDAMENTAL MODE),NOT(HIGHER MODE)')
 1020 FORMAT(' VP1,VS1,DN1 =',3F10.4)
 1030 FORMAT(' VP2,VS2,DN2 =',3F10.4)
 1040 FORMAT(' DM(DEPTH)  =',F10.4,'  DZ, DC =',2F10.4)
 1050 FORMAT(' DO  =',F10.4,'  NO, NDO =',2I10)
 1100 FORMAT(' CR2,CR1,CHECK(0) =',3D12.4)
 1200 FORMAT(//' 2*( OM, CC(I), UU(I), U/W(I), AR(I)/K )')
 1210 FORMAT(2(F8.3,2F10.4,2D11.3))
 2000 FORMAT(/' CPU TIME (SEC) =',F12.5)
 7010 FORMAT(4F10.4)
 7020 FORMAT(F10.4,3I10)
 7030 FORMAT(2(2F8.2,2D10.3))
C
      STOP
      END
C
C
C *** CAL. OF RAYLEIGH WAVE VELOCITY IN ONE LAYER MEDIUM ***
C
      REAL FUNCTION FNCCR*8(VP,VS,DZ)
      IMPLICIT REAL*8(A-H,O-Z)
C
      SP=VS/VP
      IF(SP.GE.0.0D0) Z= 0.91263D0
      IF(SP.GE.0.2D0) Z= 0.90767D0
      IF(SP.GE.0.4D0) Z= 0.88899D0
      IF(SP.GE.0.6D0) Z= 0.83575D0
      IF(SP.GE.0.7D0) Z= 0.77062D0
      IF(SP.GE.0.8D0) Z= 0.63338D0
C
      F0=FNCFR(SP,Z)
  100 Z=Z-DZ
      F1=FNCFR(SP,Z)
      AF=F0/F1
      F00=F0
      F0 =F1
      IF(AF.GE.0.D0)GOTO 100
      FZ=Z+DZ*F0/(F0-F00)
C
      FNCCR=VS*DSQRT(FZ)
C
      RETURN
      END
C
C
      REAL FUNCTION FNCFR*8(SP,Z)
      IMPLICIT REAL*8(A-H,O-Z)
C
      Z2=Z*Z
      Z3=Z2*Z
C
      F1=16.D0*(Z-1.D0)*SP*SP
      FNCFR=F1-Z3+8.D0*Z2-24.D0*Z+16.D0
C
      RETURN
      END
C
C
C *** CAL. OF EIGENVALUE EQUATION OF RAYLEIGH WAVE IN TWO LAYERS MEDIUM
C
      REAL FUNCTION FNCFR2*8(OM,C,VP1,VS1,DN1,VP2,VS2,DN2,DM)
      IMPLICIT REAL*8(A-H,O-Z)
C
      WN=OM/C
      G1=2.D0*VS1*VS1/(C*C)
      G2=2.D0*VS2*VS2/(C*C)
      CP1=C*C/(VP1*VP1)
      CP2=C*C/(VP2*VP2)
      CS1=C*C/(VS1*VS1)
      CS2=C*C/(VS2*VS2)
      P2=DSQRT(1.D0-CP2)
      S2=DSQRT(1.D0-CS2)
      P2S2=-P2*S2
C
      IF(C.LE.VP1)THEN
      X1=WN*DM*DSQRT(1.D0-CP1)
      EX1=DEXP(X1)
      EX2=DEXP(-X1)
      COP1=(EX1+EX2)/2.D0
      IF(X1.LT.1.D-5)SPP1=1.D0
      IF(X1.GE.1.D-5)SPP1=(EX1-EX2)/(2.D0*X1)
      ELSE
      P1=WN*DM*DSQRT(CP1-1.D0)
      COP1=DCOS(P1)
      IF(P1.LT.1.D-5)SPP1=1.D0
      IF(P1.GE.1.D-5)SPP1=DSIN(P1)/P1
      END IF
      IF(C.LE.VS1)THEN
      X2=WN*DM*DSQRT(1.D0-CS1)
      EX1=DEXP(X2)
      EX2=DEXP(-X2)
      COQ1=(EX1+EX2)/2.D0
      IF(X2.LT.1.D-5)SQQ1=1.0
      IF(X2.GE.1.D-5)SQQ1=(EX1-EX2)/(2.D0*X2)
      ELSE
      Q1=WN*DM*DSQRT(CS1-1.D0)
      COQ1=DCOS(Q1)
      IF(Q1.LT.1.D-5)SQQ1=1.D0
      IF(Q1.GE.1.D-5)SQQ1=DSIN(Q1)/Q1
      END IF
C
      G11=G1-1.D0
      G12=G11*G11
      G13=G12*G11
      G14=G13*G11
      G21=G2-1.D0
      G22=G21*G21
      GG=G1*G1
      G3=GG*G1
      G4=G3*G1
      P12=CP1-1.D0
      S12=CS1-1.D0
      D1=DN1/DN2
      D2=D1*D1
C
      A=G22+G2*G2*P2S2
      B=G21+G2*P2S2
      D=D1*(2.D0*G1-1.D0)*B
      E=D2*(1.D0+P2S2)*G1*G11
      C1=2.D0*G1*G11*(A-D+E)
      C2=A+C1
      C3=A*(G12+GG*P12*S12)-2.D0*B*D1*(G13+G3*P12*S12)
      C3=C3+D2*(1+P2S2)*(G14+G4*P12*S12)
      C4=D1*(G12*S2+GG*P2*S12)
      C5=D1*(G12*P2+GG*P12*S2)
C
      F2=C1-C2*COP1*COQ1+C3*WN*WN*DM*DM*SPP1*SQQ1
      FNCFR2=F2-WN*DM*(C4*COP1*SQQ1+C5*SPP1*COQ1)
C
      RETURN
      END
C
C
C *** CAL. OF GROUP VELOCITY, AMPLITUDE RATIO OF U0 & W0
C     AMPLITUDE SPECTRUM OF RAYLEIGH WAVE IN TWO LAYER MEDIUM ***
C
      SUBROUTINE UWU(OM,C,VP1,VS1,DN1,VP2,VS2,DN2,DM,U,RKL,ARK)
      IMPLICIT REAL*8(A-H,O-Z)
C
      WN=OM/C
      WD=WN*DM
      C2=C*C
      VS12=VS1*VS1
      VS22=VS2*VS2
      VP12=VP1*VP1
      VP22=VP2*VP2
      G1=2.D0*VS12/C2
      G2=2.D0*VS22/C2
      CP12=C2/VP12
      CP22=C2/VP22
      CS12=C2/VS12
      CS22=C2/VS22
      RP2=-DSQRT(1.D0-CP22)
      RS2=-DSQRT(1.D0-CS22)
C
      IF(C.LE.VP1)THEN
        P11=DSQRT(1.D0-CP12)
        X1=WD*P11
        EX1=DEXP(X1)
        EX2=DEXP(-X1)
        CP1=(EX1+EX2)/2.D0
        SP2=-P11*(EX1-EX2)/2.D0
        IF(P11.LT.1.D-5)SP1=WD
        IF(P11.GE.1.D-5)SP1=(EX1-EX2)/(2.D0*P11)
      ELSE
        P11=DSQRT(CP12-1.D0)
        P1=WD*P11
        CP1=DCOS(P1)
        SP2=P11*DSIN(P1)
        IF(P11.LT.1.D-5)SP1=WD
        IF(P11.GE.1.D-5)SP1=DSIN(P1)/P11
      END IF
C
      IF(C.LE.VS1)THEN
        Q11=DSQRT(1.D0-CS12)
        X2=WD*Q11
        EX1=DEXP(X2)
        EX2=DEXP(-X2)
        CQ1=(EX1+EX2)/2.D0
        SQ2=-Q11*(EX1-EX2)/2.D0
        IF(Q11.LT.1.D-5)SQ1=WD
        IF(Q11.GE.1.D-5)SQ1=(EX1-EX2)/(2.D0*Q11)
      ELSE
        Q11=DSQRT(CS12-1.D0)
        Q1=WD*Q11
        CQ1=DCOS(Q1)
        SQ2=Q11*DSIN(Q1)
        IF(Q11.LT.1.D-5)SQ1=WD
        IF(Q11.GE.1.D-5)SQ1=DSIN(Q1)/Q11
      END IF
C
      D12=DN1/DN2
      G11=G1
      G12=G1-1.D0
      G13=G1*G1
      G14=G12*G12
      G15=DN1*C2
      G16=DN2*C2
C
      C11=RP2
      C12=RS2
      C13=G2-1.D0
      C14=G2*C11
      C15=G2*C12
C
      A11=G11*CP1-G12*CQ1
      A12=G12*SP1+G11*SQ2
      A13=-(CP1-CQ1)/G15
      A14=(SP1+SQ2)/G15
      A21=-G11*SP2-G12*SQ1
      A22=-G12*CP1+G11*CQ1
      A23=(SP2+SQ1)/G15
      A24=A13
      A31=D12*G11*G12*(CP1-CQ1)
      A32=D12*(G14*SP1+G13*SQ2)
      A33=A22/G16
      A34=A12/G16
      A41=D12*(G13*SP2+G14*SQ1)
      A42=A31
      A43=A21/G16
      A44=A11/G16
C
      RL= C14*A11+C13*A21-C11*A31+A41
      RK=-C14*A12+C13*A22+C11*A32+A42
      RG= C14*A13+C13*A23-C11*A33+A43
      RR=-C14*A14+C13*A24+C11*A34+A44
      RN=-C13*A11-C15*A21+A31-C12*A41
      RM=-C13*A12+C15*A22+A32+C12*A42
      RH=-C13*A13-C15*A23+A33-C12*A43
      RS=-C13*A14+C15*A24+A34+C12*A44
C
      RP12=CP12-1.D0
      RS12=CS12-1.D0
      CD=C*DM
C
      G1K=2.D0*G1/WN
      G2K=2.D0*G2/WN
      CK=-C/WN
      RP2K=-C*CK/VP22/RP2
      RS2K=-C*CK/VS22/RS2
      RP11K=-DM/RP12
      RS11K=-DM/RS12
      RP12K=-DM
      RS12K=-DM
      SP11K=-DM*SP1
      SS11K=-DM*SQ1
      SP12K=-CP12*SP1/WN
      SS12K=-CS12*SQ1/WN
C
      G11K=G1K
      G12K=G1K
      G13K=2.D0*G1*G1K
      G14K=2.D0*G12*G1K
C
      CP1K=-SP11K
      CQ1K=-SS11K
      SP1K=CP1*RP11K-SP12K/RP12
      SQ1K=CQ1*RS11K-SS12K/RS12
      SP2K=SP12K+CP1*RP12K
      SQ2K=SS12K+CQ1*RS12K
C
      C11K=RP2K
      C12K=RS2K
      C13K=G2K
      C14K=G2K*C11+G2*C11K
      C15K=G2K*C12+G2*C12K
C
      A11K= G11K*CP1+G11*CP1K-G12K*CQ1-G12*CQ1K
      A12K= G12K*SP1+G12*SP1K+G11K*SQ2+G11*SQ2K
      A21K=-G11K*SP2-G11*SP2K-G12K*SQ1-G12*SQ1K
      A22K=-G12K*CP1-G12*CP1K+G11K*CQ1+G11*CQ1K
      A31K=D12*((G11K*G12+G11*G12K)*(CP1-CQ1)+G11*G12*(CP1K-CQ1K))
      A32K=D12*(G14K*SP1+G14*SP1K+G13K*SQ2+G13*SQ2K)
      A41K=D12*(G13K*SP2+G13*SP2K+G14K*SQ1+G14*SQ1K)
      A42K=A31K
C
      RLK= C14K*A11+C14*A11K+C13K*A21+C13*A21K
      RLK=RLK-C11K*A31-C11*A31K+A41K
      RKK=-C14K*A12-C14*A12K+C13K*A22+C13*A22K
      RKK=RKK+C11K*A32+C11*A32K+A42K
      RNK=-C13K*A11-C13*A11K-C15K*A21-C15*A21K
      RNK=RNK+A31K-C12K*A41-C12*A41K
      RMK=-C13K*A12-C13*A12K+C15K*A22+C15*A22K
      RMK=RMK+A32K+C12K*A42+C12*A42K
C
      FRK=RNK*RK+RN*RKK+RLK*RM+RL*RMK
C
      G1W=-2.D0*G1/WN/C
      G2W=-2.D0*G2/WN/C
      CW=1.D0/WN
      RP2W=-C*CW/VP22/RP2
      RS2W=-C*CW/VS22/RS2
      RP11W=CD/VP12/RP12
      RS11W=CD/VS12/RS12
      RP12W=CD/VP12
      RS12W=CD/VS12
      SP11W=RP12W*SP1
      SS11W=RS12W*SQ1
      SP12W=C*SP1/VP12/WN
      SS12W=C*SQ1/VS12/WN
C
      G11W=G1W
      G12W=G1W
      G13W=2.D0*G1*G1W
      G14W=2.D0*G12*G1W
C
      CP1W=-SP11W
      CQ1W=-SS11W
      SP1W=CP1*RP11W-SP12W/RP12
      SQ1W=CQ1*RS11W-SS12W/RS12
      SP2W=SP12W+CP1*RP12W
      SQ2W=SS12W+CQ1*RS12W
C
      C11W=RP2W
      C12W=RS2W
      C13W=G2W
      C14W=G2W*C11+G2*C11W
      C15W=G2W*C12+G2*C12W
C
      A11W= G11W*CP1+G11*CP1W-G12W*CQ1-G12*CQ1W
      A12W= G12W*SP1+G12*SP1W+G11W*SQ2+G11*SQ2W
      A21W=-G11W*SP2-G11*SP2W-G12W*SQ1-G12*SQ1W
      A22W=-G12W*CP1-G12*CP1W+G11W*CQ1+G11*CQ1W
      A31W=D12*((G11W*G12+G11*G12W)*(CP1-CQ1)+G11*G12*(CP1W-CQ1W))
      A32W=D12*(G14W*SP1+G14*SP1W+G13W*SQ2+G13*SQ2W)
      A41W=D12*(G13W*SP2+G13*SP2W+G14W*SQ1+G14*SQ1W)
      A42W=A31W
C
      RLW= C14W*A11+C14*A11W+C13W*A21+C13*A21W
      RLW=RLW-C11W*A31-C11*A31W+A41W
      RKW=-C14W*A12-C14*A12W+C13W*A22+C13*A22W
      RKW=RKW+C11W*A32+C11*A32W+A42W
      RNW=-C13W*A11-C13*A11W-C15W*A21-C15*A21W
      RNW=RNW+A31W-C12W*A41-C12*A41W
      RMW=-C13W*A12-C13*A12W+C15W*A22+C15*A22W
      RMW=RMW+A32W+C12W*A42+C12*A42W
C
      FRW=RNW*RK+RN*RKW+RLW*RM+RL*RMW
C
      U=-FRK/FRW
      RKL=RK/RL
      RMN=RM/RN
      FR =RN*RK+RL*RM
      ARK=(RG*RN-RL*RH)/FRK/WN
      RGSH=(RR*RN+RS*RL)/(RG*RN-RH*RL)
      CHECK1=DABS(-RKL/RMN-1.D0)
      CHECK2=DABS(RGSH/RKL-1.D0)
      IF(CHECK1.GT.1.D-3)THEN
      WRITE(*,*)'OM,CHECK1 !!!!! FROM SUB.UWU',OM,CHECK1
      WRITE(*,*)'OM,CHECK2 !!!!! FROM SUB.UWU',OM,CHECK2
      END IF
C
      RETURN
      END
