C
C *** GROUP VELOCITY AND AMPLITUDE SPECTRA OF RAYLEIGH WAVES 
C     DUE TO HORIZONTAL AND VERTICAL POINT FORCE IN A MULTI-LAYERED MEDIUM ***
C
      SUBROUTINE RGRAMP(DOM,NOM,MOM,ML,NL,DNS,VP,VS,DEP
     *                 ,CC,NMODE,CFOM,NSMD,MMD
     *                 ,NLS,SDEP,UU,UW0,ARV,ARH)
C
C ** INPUT DATA **
C DOM     : DELTA OMEGA
C NOM     : TOTAL NUMBER OF OMEGA
C MOM     : NOM IN MAIN PROGRAM DIMENSION
C ML      : NL IN MAIN
C NL      : SOIL LAYERS NUMBER
C DNS(I)  : DENSITY IN Ith LAYER
C VP(I)   : VP IN Ith LAYER
C VS(I)   : VS IN Ith LAYER
C DEP(I)  : THICKNESS OF Ith LAYER
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 NLS     : SOURCE LOCATED LAYER
C SDEP    : SOURCE DEPTH FROM THE UPPER BOUNDARY OF THE SOURCE LAYER
C ** OUTPUT DATA **
C UU(I,J) : GROUP VELOCITY OF Jth MODE AT Ith OMEGA
C UW0(I,J): AMPLITUDE RATIO ( IMAG(U0/W0) ) ON SURFACE OF Jth MODE 
C           AT Ith OMEGA
C ARV(I,J): AMPLITUDE DISPLACEMENT SPECTRUM DUE TO VERTICAL UNIT POINT SOURCE
C ARH(I,J): AMPLITUDE DISPLACEMENT SPECTRUM DUE TO 
C           HORIZONTAL UNIT POINT SOURCE (IMAG(ARH))
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION VP(ML),VS(ML),DEP(ML),DNS(ML)
     *         ,CC(MOM,0:MMD),NMODE(MOM),UU(MOM,0:MMD)
     *         ,CFOM(0:MOM),UW0(MOM,0:MMD)
     *         ,ARV(MOM,0:MMD),ARH(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
     *          ,SB13,SB23,SB33,SB43,SB14,SB24,SB34,SB44
     *          ,SRL,SRK,SRN,SRM,SRG,SRR,SRH,SRS
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=0,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       * 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)
          B13=-RAN*B14
C
          B21=-B12
          B22= GN*RBN
          B23= B14
          B24= RBN*B14
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 IL=NL-1,1,-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 & SDQ AND CP, STP &SDP *
C
            BM=VS(IL)
            BM2=BM*BM
            CB2=CV*CV/BM2
            CALL RYSC(WN,BM,CV,CB2,DEP(IL),CQ,STQ,SDQ)
C
            AM=VP(IL)
            AM2=AM*AM
            CA2=CV*CV/AM2
            CALL RYSC(WN,AM,CV,CA2,DEP(IL),CP,STP,SDP)
C
C         * CPK,CPW,-,SDQK,SDQW *
C
            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(IL)*(STP+WN*CVAK*SDP)
            CPW=-DEP(IL)*WN*CVAW*SDP
            CQK=-DEP(IL)*(STQ+WN*CVBK*SDQ)
            CQW=-DEP(IL)*WN*CVBW*SDQ
            STPK=CVAK*SDP+DEP(IL)*CP*(RA2+WN*CVAK)
            STPW=CVAW*(SDP+DEP(IL)*CP*WN)
            STQK=CVBK*SDQ+DEP(IL)*CQ*(RB2+WN*CVBK)
            STQW=CVBW*(SDQ+DEP(IL)*CQ*WN)
            SDPK=-CVAK*SDP/RA2+DEP(IL)*CP*(1.D0+WN*CVAK/RA2)
            SDPW=CVAW*(-SDP+DEP(IL)*CP*WN)/RA2
            SDQK=-CVBK*SDQ/RB2+DEP(IL)*CQ*(1.D0+WN*CVBK/RB2)
            SDQW=CVBW*(-SDQ+DEP(IL)*CQ*WN)/RB2
C
C         * C11-C16, D11-D15 AND SB MATRIX *
C
            GM=2.D0/CB2
            CALL RYSB(GM,DNS(IL),CV,CP,CQ,SDP,SDQ,STP,STQ
     *               ,C11,C12,C13,C14,C15,C16,D11,D12,D13,D14,D15,SB)
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(IL)*CV*CK
            C16W=2.D0*DNS(IL)*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         * SGS - SSS
C
            IF(IL.EQ.NLS)THEN
C
              CALL RYSC(WN,BM,CV,CB2,DEP(IL)-SDEP,SCQ,SSTQ,SSDQ)
              CALL RYSC(WN,AM,CV,CA2,DEP(IL)-SDEP,SCP,SSTP,SSDP)
C
              SD11=SCP-SCQ
              SD12=SSDP+SSTQ
              SD13=SSTP+SSDQ
C
              SB13=-DCMPLX(C15*SD11,0.D0)
              SB23= DCMPLX(0.D0,C15*SD13)
              SB33= DCMPLX(-C12*SCP+C11*SCQ,0.D0)
              SB43=-DCMPLX(0.0D0,C11*SSTP+C12*SSDQ)
              SB14= DCMPLX(0.D0,C15*SD12)
              SB24= SB13
              SB34= DCMPLX(0.0D0,C12*SSDP+C11*SSTQ)
              SB44= DCMPLX(C11*SCP-C12*SCQ,0.0D0)
C
              SRL=B11*AA(1,1)+B12*AA(2,1)+B13*AA(3,1)+B14*AA(4,1)
              SRK=B11*AA(1,2)+B12*AA(2,2)+B13*AA(3,2)+B14*AA(4,2)
              SRG=B11*AA(1,3)+B12*AA(2,3)+B13*AA(3,3)+B14*AA(4,3)
              SRR=B11*AA(1,4)+B12*AA(2,4)+B13*AA(3,4)+B14*AA(4,4)
C
              SRN=B21*AA(1,1)+B22*AA(2,1)+B23*AA(3,1)+B24*AA(4,1)
              SRM=B21*AA(1,2)+B22*AA(2,2)+B23*AA(3,2)+B24*AA(4,2)
              SRH=B21*AA(1,3)+B22*AA(2,3)+B23*AA(3,3)+B24*AA(4,3)
              SRS=B21*AA(1,4)+B22*AA(2,4)+B23*AA(3,4)+B24*AA(4,4)
C
              GSS=IMAG(SRL*SB13+SRK*SB23+SRG*SB33+SRR*SB43)
              HSS=REAL(SRN*SB13+SRM*SB23+SRH*SB33+SRS*SB43)
              RSS=REAL(SRL*SB14+SRK*SB24+SRG*SB34+SRR*SB44)
              SSS=IMAG(SRN*SB14+SRM*SB24+SRH*SB34+SRS*SB44)
C
            END IF
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) +SA(I,K)*SB(K,J)
                  AAK(I,J)=AAK(I,J)+SAK(I,K)*SB(K,J)+SA(I,K)*SBK(K,J)
                  AAW(I,J)=AAW(I,J)+SAW(I,K)*SB(K,J)+SA(I,K)*SBW(K,J)
  270           CONTINUE
  260         CONTINUE
  250       CONTINUE
C
  200     CONTINUE
C
C      ** RL - RM *
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 ,UW0 AND ARV & ARH **
C
          UU(IOM,MODE) =-FRK/FRW
C         ARK=(RGI*RNR-RLI*RHR)/(WN*FRK)
          UW0(IOM,MODE)=-IMAG(RK/RL)
C
          IF(NLS.EQ.NL)THEN
            RRA=-IMAG(RAN)
            RRB=-IMAG(RBN)
            EXA=DEXP(-WN*RRA*SDEP)
            EXB=DEXP(-WN*RRB*SDEP)
            GSS= B14*RRA*EXA
            HSS= B14*EXB
            RSS= B14*EXA
            SSS=-B14*RRB*EXB
          END IF
C
          SB1=GSS*RNR-HSS*RLI
          SB2=RSS*RNR+SSS*RLI
C
          ARV(IOM,MODE)=-SB1/(WN*FRK)
          ARH(IOM,MODE)= SB2/(WN*FRK)
            
  110   CONTINUE
C
  100 CONTINUE   
C
      RETURN
      END 
