C
C *** RAYLEIGH WAVE GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
      SUBROUTINE RAYUA(DOM,NOM,MOM,CC,NMODE,CFOM,NSMD,MMD,UU,ARK,UW0)
C
C ** INPUT DATA **
C DOM     : DELTA OMEGA
C NOM     : TOTAL NUMBER OF OMEGA
C MOM     : NOM IN MAIN PROGRAM DIMENSION
C CC(I,J) : PHASE VELOCITY OF Jth MODE AT Ith OMEGA
C NMODE(I): OUTPUT HIGHEST MODE AT Ith OMEGA
C CFOM(I) : CUT-OFF OMEGA OF Ith MODE
C NSMD    : START MODE
C MMD     : NEMD IN MAIN PROGRAM DIMENSION
C ** OUTPUT DATA **
C UU(I,J) : GROUP VELOCITY OF Jth MODE AT Ith OMEGA
C ARK(I,J): AMPLITUDE SPECTRUM (AR/K) OF Jth RAYLEIGH MODE AT Ith OMEGA
C UW0(I,J): AMPLITUDE RATIO ( IMAG(U0/W0) ) ON SURFACE OF Jth MODE 
C           AT Ith OMEGA
C
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON / LAYER / NL,VP(40),VS(40),DEP(40),SM(40)
     *       / DNSITY / DNS(40)
      DIMENSION CC(MOM,0:MMD),NMODE(MOM),UU(MOM,0:MMD),ARK(MOM,0:MMD)
     *         ,CFOM(0:MOM),UW0(MOM,0:MMD)
      COMPLEX*16 AA(4,4),SA(4,4),SB(4,4),AAK(4,4),SAK(4,4),SBK(4,4)
     *          ,AAW(4,4),SAW(4,4),SBW(4,4)
     *          ,RAN,RBN,B11,B22,B13,B24,RL,RK,RN,RM
     *          ,RANK,RBNK,B11K,B22K,B13K,B24K,RLK,RKK,RNK,RMK
     *          ,RANW,RBNW,B11W,B22W,B13W,B24W,RLW,RKW,RNW,RMW
C
C *** OMEGA ***
C
      DO 100 IOM=1,NOM
C
        OM=DOM*DFLOAT(IOM)
        IF(OM.LT.CFOM(NSMD))GOTO 100
C
C *** MODE ***
C
        DO 110 MODE=NSMD,NMODE(IOM)

C      ** WAVE NUMBER **
C
          CV=CC(IOM,MODE)
          WN=OM/CV
C
C      ** CK AND CW **
C
          CK=-CV/WN
          CW=1.D0/WN
C
C      ** INITIALIZATION OF AA,AAK AND AAW MATRIX **
C
          DO 120 I=1,4
            DO 130 J=1,4
              AA(I,J) =(0.0D0,0.0D0)
              AAK(I,J)=(0.0D0,0.0D0)
              AAW(I,J)=(0.0D0,0.0D0)
  130       CONTINUE
            AA(I,I)=(1.0D0,0.0D0)
  120     CONTINUE
C
C      ** AA MATRIX **
C
          DO 200 ML=1,NL-1
C
            DO 210 I=1,4
              DO 220 J=1,4
                SA(I,J) =AA(I,J)
                SAK(I,J)=AAK(I,J)
                SAW(I,J)=AAW(I,J)
  220         CONTINUE
  210       CONTINUE
C
C         * CQ, STQ AND SDQ *
C
            BM=VS(ML)
            BM2=BM*BM
            CB2=CV*CV/BM2
            IF(BM-CV)10,20,30
C
   10         RBM=DSQRT(CB2-1.D0)
              QM=WN*RBM*DEP(ML)
              CQ=DCOS(QM)
              SQ=DSIN(QM)
              STQ=RBM*SQ
              SDQ=SQ/RBM
              GOTO 40
C
   20         WRITE(*,*)'FROM SUBROUTINE RAYUA :'
              WRITE(*,*)'CV=BM THEN STOP !!!'
              WRITE(*,*)'OM,MODE,CV,BM =',OM,MODE,CV,BM
              STOP
C
   30         RBM=DSQRT(1.D0-CB2)
              QM=WN*RBM*DEP(ML)
              EX1=DEXP(QM)
              EX2=DEXP(-QM)
              CQ=(EX1+EX2)/2.0D0
              EX12=(EX1-EX2)/2.0D0
              STQ=-RBM*EX12
              SDQ=EX12/RBM
C
C         * CP, STP AND SDP *
C
   40       AM=VP(ML)
            AM2=AM*AM
            CA2=CV*CV/AM2
            IF(AM-CV)50,60,70
C
   50         RAM=DSQRT(CA2-1.D0)
              PM=WN*RAM*DEP(ML)
              CP=DCOS(PM)
              SP=DSIN(PM)
              STP=RAM*SP
              SDP=SP/RAM
              GOTO 80
C
   60         WRITE(*,*)'FROM SUBROUTINE RAYUA :'
              WRITE(*,*)'CV=AM THEN STOP !!!'
              WRITE(*,*)'OM,MODE,CV,AM =',OM,MODE,CV,AM
              STOP
C
   70         RAM=DSQRT(1.D0-CA2)
              PM=WN*RAM*DEP(ML)
              EX1=DEXP(PM)
              EX2=DEXP(-PM)
              CP=(EX1+EX2)/2.0D0
              EX12=(EX1-EX2)/2.0D0
              STP=-RAM*EX12
              SDP=EX12/RAM
C
C         * CPK,CPW,-,SDQK,SDQW *
C
   80       CVAK=CV*CK/AM2
            CVAW=CV*CW/AM2
            CVBK=CV*CK/BM2
            CVBW=CV*CW/BM2
C
            RB2=CB2-1.D0
            RA2=CA2-1.D0
            CPK=-DEP(ML)*(STP+WN*CVAK*SDP)
            CPW=-DEP(ML)*WN*CVAW*SDP
            CQK=-DEP(ML)*(STQ+WN*CVBK*SDQ)
            CQW=-DEP(ML)*WN*CVBW*SDQ
            STPK=CVAK*SDP+DEP(ML)*CP*(RA2+WN*CVAK)
            STPW=CVAW*(SDP+DEP(ML)*CP*WN)
            STQK=CVBK*SDQ+DEP(ML)*CQ*(RB2+WN*CVBK)
            STQW=CVBW*(SDQ+DEP(ML)*CQ*WN)
            SDPK=-CVAK*SDP/RA2+DEP(ML)*CP*(1.D0+WN*CVAK/RA2)
            SDPW=CVAW*(-SDP+DEP(ML)*CP*WN)/RA2
            SDQK=-CVBK*SDQ/RB2+DEP(ML)*CQ*(1.D0+WN*CVBK/RB2)
            SDQW=CVBW*(-SDQ+DEP(ML)*CQ*WN)/RB2
C
C         * C11-C16, D11-D15 AND SB MATRIX *
C
            GM=2.D0/CB2
C
            C11=GM
            C12=GM-1.D0
            C13=GM*GM
            C14=C12*C12
            C16=DNS(ML)*CV*CV
            C15=1.D0/C16
C
            D11=CP-CQ
            D12=SDP+STQ
            D13=STP+SDQ
            D14=C14*SDP+C13*STQ
            D15=C13*STP+C14*SDQ
C
            SB(1,1)= DCMPLX(C11*CP-C12*CQ,0.0D0)
            SB(1,2)= DCMPLX(0.0D0,C12*SDP+C11*STQ)
            SB(1,3)=-DCMPLX(C15*D11,0.D0)
            SB(1,4)= DCMPLX(0.D0,C15*D12)
            SB(2,1)=-DCMPLX(0.0D0,C11*STP+C12*SDQ)
            SB(2,2)= DCMPLX(-C12*CP+C11*CQ,0.D0)
            SB(2,3)= DCMPLX(0.D0,C15*D13)
            SB(2,4)= SB(1,3)
            SB(3,1)= DCMPLX(C16*C11*C12*D11,0.0D0)
            SB(3,2)= DCMPLX(0.D0,C16*D14)
            SB(3,3)= SB(2,2)
            SB(3,4)= SB(1,2)
            SB(4,1)= DCMPLX(0.D0,C16*D15)
            SB(4,2)= SB(3,1)
            SB(4,3)= SB(2,1)
            SB(4,4)= SB(1,1)
C
C         * C11K-C16W, D11K-D15W, SBK AND SBW MATRIX *
C
            GMK=-2.D0*GM*CK/CV
            GMW=-2.D0*GM*CW/CV
C
            C11K=GMK
            C11W=GMW
            C12K=GMK
            C12W=GMW
            C13K=2.D0*GM*GMK
            C13W=2.D0*GM*GMW
            C14K=2.D0*C12*C12K
            C14W=2.D0*C12*C12W
            C16K=2.D0*DNS(ML)*CV*CK
            C16W=2.D0*DNS(ML)*CV*CW
            C15K=-C16K/(C16*C16)
            C15W=-C16W/(C16*C16)
C
            D11K=CPK-CQK
            D11W=CPW-CQW
            D12K=SDPK+STQK
            D12W=SDPW+STQW
            D13K=STPK+SDQK
            D13W=STPW+SDQW
            D14K=C14K*SDP+C14*SDPK+C13K*STQ+C13*STQK
            D14W=C14W*SDP+C14*SDPW+C13W*STQ+C13*STQW
            D15K=C13K*STP+C13*STPK+C14K*SDQ+C14*SDQK
            D15W=C13W*STP+C13*STPW+C14W*SDQ+C14*SDQW
C
            SBK(1,1)= DCMPLX(C11K*CP+C11*CPK-C12K*CQ-C12*CQK,0.0D0)
            SBK(1,2)= DCMPLX(0.0D0,C12K*SDP+C12*SDPK+C11K*STQ+C11*STQK)
            SBK(1,3)=-DCMPLX(C15K*D11+C15*D11K,0.D0)
            SBK(1,4)= DCMPLX(0.D0,C15K*D12+C15*D12K)
            SBK(2,1)=-DCMPLX(0.0D0,C11K*STP+C11*STPK+C12K*SDQ+C12*SDQK)
            SBK(2,2)= DCMPLX(-C12K*CP-C12*CPK+C11K*CQ+C11*CQK,0.D0)
            SBK(2,3)= DCMPLX(0.D0,C15K*D13+C15*D13K)
            SBK(2,4)= SBK(1,3)
            SBK(3,1)= DCMPLX((C16K*C11+C16*C11K)*C12*D11,0.0D0)
            SBK(3,1)= SBK(3,1)+DCMPLX(C16*C11*(C12K*D11+C12*D11K),0.0D0)
            SBK(3,2)= DCMPLX(0.D0,C16K*D14+C16*D14K)
            SBK(3,3)= SBK(2,2)
            SBK(3,4)= SBK(1,2)
            SBK(4,1)= DCMPLX(0.D0,C16K*D15+C16*D15K)
            SBK(4,2)= SBK(3,1)
            SBK(4,3)= SBK(2,1)
            SBK(4,4)= SBK(1,1)
C
            SBW(1,1)= DCMPLX(C11W*CP+C11*CPW-C12W*CQ-C12*CQW,0.0D0)
            SBW(1,2)= DCMPLX(0.0D0,C12W*SDP+C12*SDPW+C11W*STQ+C11*STQW)
            SBW(1,3)=-DCMPLX(C15W*D11+C15*D11W,0.D0)
            SBW(1,4)= DCMPLX(0.D0,C15W*D12+C15*D12W)
            SBW(2,1)=-DCMPLX(0.0D0,C11W*STP+C11*STPW+C12W*SDQ+C12*SDQW)
            SBW(2,2)= DCMPLX(-C12W*CP-C12*CPW+C11W*CQ+C11*CQW,0.D0)
            SBW(2,3)= DCMPLX(0.D0,C15W*D13+C15*D13W)
            SBW(2,4)= SBW(1,3)
            SBW(3,1)= DCMPLX((C16W*C11+C16*C11W)*C12*D11,0.0D0)
            SBW(3,1)= SBW(3,1)+DCMPLX(C16*C11*(C12W*D11+C12*D11W),0.0D0)
            SBW(3,2)= DCMPLX(0.D0,C16W*D14+C16*D14W)
            SBW(3,3)= SBW(2,2)
            SBW(3,4)= SBW(1,2)
            SBW(4,1)= DCMPLX(0.D0,C16W*D15+C16*D15W)
            SBW(4,2)= SBW(3,1)
            SBW(4,3)= SBW(2,1)
            SBW(4,4)= SBW(1,1)
            
C
C         * AA MATRIX *
C
            DO 250 I=1,4
              DO 260 J=1,4
                AA(I,J)=(0.D0,0.D0)
                AAK(I,J)=(0.D0,0.D0)
                AAW(I,J)=(0.D0,0.D0)
                DO 270 K=1,4
                  AA(I,J) =AA(I,J) +SB(I,K)*SA(K,J)
                  AAK(I,J)=AAK(I,J)+SBK(I,K)*SA(K,J)+SB(I,K)*SAK(K,J)
                  AAW(I,J)=AAW(I,J)+SBW(I,K)*SA(K,J)+SB(I,K)*SAW(K,J)
  270           CONTINUE
  260         CONTINUE
  250       CONTINUE
C
  200     CONTINUE
C
C       * GN,RAN,RBN,B11-B24 AND RL-RM ***
C
          GN=2.D0*VS(NL)*VS(NL)/(CV*CV)
          RAN=-DCMPLX(0.D0,DSQRT(1.0D0-CV*CV/(VP(NL)*VP(NL))))
          RBN=-DCMPLX(0.D0,DSQRT(1.0D0-CV*CV/(VS(NL)*VS(NL))))
C
          B11=GN*RAN
          B12=GN-1.D0
          B14=1.D0/(DNS(NL)*CV*CV)
          B22=GN*RBN
          B13=-RAN*B14
          B21=-B12
          B23=B14
          B24=RBN*B14
C
          RL=B11*AA(1,1)+B12*AA(2,1)+B13*AA(3,1)+B14*AA(4,1)
          RK=B11*AA(1,2)+B12*AA(2,2)+B13*AA(3,2)+B14*AA(4,2)
          RN=B21*AA(1,1)+B22*AA(2,1)+B23*AA(3,1)+B24*AA(4,1)
          RM=B21*AA(1,2)+B22*AA(2,2)+B23*AA(3,2)+B24*AA(4,2)
C
C       * GNK,RANK,RBNK,B11K-B24K AND RLK-RMK ***
C
          GNK=-2.D0*GN*CK/CV
          RANK=CV*CK/(VP(NL)*VP(NL)*RAN)
          RBNK=CV*CK/(VS(NL)*VS(NL)*RBN)
C
          B11K=GNK*RAN+GN*RANK
          B12K=GNK
          B14K=-2.D0*CK*B14/CV
          B22K=GNK*RBN+GN*RBNK
          B13K=-RANK*B14-RAN*B14K
          B21K=-B12K
          B23K=B14K
          B24K=RBNK*B14+RBN*B14K
C
          RLK=B11K*AA(1,1)+B11*AAK(1,1)+B12K*AA(2,1)+B12*AAK(2,1)
          RLK=RLK+B13K*AA(3,1)+B13*AAK(3,1)+B14K*AA(4,1)+B14*AAK(4,1)
          RKK=B11K*AA(1,2)+B11*AAK(1,2)+B12K*AA(2,2)+B12*AAK(2,2)
          RKK=RKK+B13K*AA(3,2)+B13*AAK(3,2)+B14K*AA(4,2)+B14*AAK(4,2)
          RNK=B21K*AA(1,1)+B21*AAK(1,1)+B22K*AA(2,1)+B22*AAK(2,1)
          RNK=RNK+B23K*AA(3,1)+B23*AAK(3,1)+B24K*AA(4,1)+B24*AAK(4,1)
          RMK=B21K*AA(1,2)+B21*AAK(1,2)+B22K*AA(2,2)+B22*AAK(2,2)
          RMK=RMK+B23K*AA(3,2)+B23*AAK(3,2)+B24K*AA(4,2)+B24*AAK(4,2)
C
C       * GNW,RANW,RBNW,B11W-B24W AND RLW-RMW ***
C
          GNW=-2.D0*GN*CW/CV
          RANW=CV*CW/(VP(NL)*VP(NL)*RAN)
          RBNW=CV*CW/(VS(NL)*VS(NL)*RBN)
C
          B11W=GNW*RAN+GN*RANW
          B12W=GNW
          B14W=-2.D0*CW*B14/CV
          B22W=GNW*RBN+GN*RBNW
          B13W=-RANW*B14-RAN*B14W
          B21W=-B12W
          B23W=B14W
          B24W=RBNW*B14+RBN*B14W
C
          RLW=B11W*AA(1,1)+B11*AAW(1,1)+B12W*AA(2,1)+B12*AAW(2,1)
          RLW=RLW+B13W*AA(3,1)+B13*AAW(3,1)+B14W*AA(4,1)+B14*AAW(4,1)
          RKW=B11W*AA(1,2)+B11*AAW(1,2)+B12W*AA(2,2)+B12*AAW(2,2)
          RKW=RKW+B13W*AA(3,2)+B13*AAW(3,2)+B14W*AA(4,2)+B14*AAW(4,2)
          RNW=B21W*AA(1,1)+B21*AAW(1,1)+B22W*AA(2,1)+B22*AAW(2,1)
          RNW=RNW+B23W*AA(3,1)+B23*AAW(3,1)+B24W*AA(4,1)+B24*AAW(4,1)
          RMW=B21W*AA(1,2)+B21*AAW(1,2)+B22W*AA(2,2)+B22*AAW(2,2)
          RMW=RMW+B23W*AA(3,2)+B23*AAW(3,2)+B24W*AA(4,2)+B24*AAW(4,2)
C
C      ** FRK,FRW, RG AND RH **
C
          FRK=REAL(RNK*RK+RN*RKK-RLK*RM-RL*RMK)
          FRW=REAL(RNW*RK+RN*RKW-RLW*RM-RL*RMW)
C
          RGI=IMAG(B11*AA(1,3)+B12*AA(2,3)+B13*AA(3,3)+B14*AA(4,3))
          RNR=REAL(RN)
          RLI=IMAG(RL)
          RHR=REAL(B21*AA(1,3)+B22*AA(2,3)+B23*AA(3,3)+B24*AA(4,3))
C
C      ** UU AND ARK **
C
          UU(IOM,MODE) =-FRK/FRW
          ARK(IOM,MODE)=(RGI*RNR-RLI*RHR)/(WN*FRK)
          UW0(IOM,MODE)=-IMAG(RK/RL)
C
  110   CONTINUE
C
  100  CONTINUE
C                     
      RETURN
      END
