C
C *** RAYLEIGH WAVE GROUP VELOCITY AND AMPLITUDE SPECTRUM ***
C
      SUBROUTINE RAYKIN(DOM,NOM,MOM,ML,NL,DNS,VP,VS,DEP
     *                 ,CC,NMODE,CFOM,NSMD,MMD,UU,ARK,UW0
     *                 ,NFSL,NFNL,DDEP,NNL,U0,W0
     *                 ,DGS,DHS,DRS,DSS,DGSZ,DHSZ,DRSZ,DSSZ)
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 NSMD    : START MODE
C MMD     : NEMD IN MAIN PROGRAM DIMENSION
C NFSL    : FIRST LAYER FOR U0 AND W0
C NFNL    : FINAL LAYER FOR U0 AND W0
C DDEP    : DELTA THICKNESS FOR Nth LAYER INTEGRATION
C NNL     : INTEGRAL NUMBER FOR DDEP
C DGS(I) - DSSZ(I) : WORK SPACE
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 U0(I)   : APPROXIMATE RAYLEIGH WAVE(HORIZONTAL DISPLACEMENT
C           ON FREE-SURFACE) IN A SEDIMENTARY BASIN 
C           DUE TO 2-D SV HASKELL'S FAULT MODEL
C W0(I)   : APPROXIMATE RAYLEIGH WAVE(VERTICAL DISPLACEMENT
C           ON FREE-SURFACE) IN A SEDIMENTARY BASIN 
C           DUE TO 2-D SV HASKELL'S FAULT MODEL
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION VP(ML),VS(ML),DEP(ML),DNS(ML)
     *         ,DGS(ML),DHS(ML),DRS(ML),DSS(ML)
     *         ,DGSZ(ML),DHSZ(ML),DRSZ(ML),DSSZ(ML)
     *         ,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,RG,RR,RH,RS
     *          ,RANK,RBNK,B11K,B22K,B13K,B24K,RLK,RKK,RNK,RMK
     *          ,RANW,RBNW,B11W,B22W,B13W,B24W,RLW,RKW,RNW,RMW
     *          ,DB13,DB23,DB33,DB43,DB14,DB24,DB34,DB44
     *          ,DB13Z,DB23Z,DB33Z,DB43Z,DB14Z,DB24Z,DB34Z,DB44Z
     *          ,U0(MOM),W0(MOM),DUP,DWP,DUS,DWS,UF,WF
     *          ,DUZZ,DUZX,DUXZ,DUXX,DWZZ,DWZX,DWXZ,DWXX
     *          ,DUGN,DUHN,DURN,DUSN,DWGN,DWHN,DWRN,DWSN
     *          ,DUGZ,DUHZ,DURZ,DUSZ,DWGZ,DWHZ,DWRZ,DWSZ
     *          ,DUB1,DUB2,DUB1Z,DUB2Z,DWB1,DWB2,DWB1Z,DWB2Z
     *          ,DUUZZ,DUUZX,DUUXZ,DUUXX,DUWZZ,DUWZX,DUWXZ,DUWXX
     *          ,DWUZZ,DWUZX,DWUXZ,DWUXX,DWWZZ,DWWZX,DWWXZ,DWWXX
C
C *** PARAMETERS INPUT FOR 2D-SV HASKELL'S FAULT MODEL ***
C
      CALL HSKIN(XF,ZF,DS0,P0,S0,WD,DIP,DIS,RV,RT) 
C
C *** OMEGA ***
C
      DO 100 IOM=1,NOM
C
        OM=DOM*DFLOAT(IOM)
        IF(OM.LT.CFOM(NSMD))GOTO 100
C
C    ** INITIALIZATION OF U0 AND W0 **
C
        U0(IOM)=(0.D0,0.D0)
        W0(IOM)=(0.D0,0.D0)
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
   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(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
C
            C11=GM
            C12=GM-1.D0
            C13=GM*GM
            C14=C12*C12
            C16=DNS(IL)*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(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         * DCP - DSDQ *
C
            DCP =-SDP/WN
            DCQ =-SDQ/WN
            DSTP= (1.D0-CP)/WN
            DSTQ= (1.D0-CQ)/WN
            DSDP= DSTP/RA2
            DSDQ= DSTQ/RB2
C
C         * DB13 - DB44 *
C
            DB13=DCMPLX(-C15*(DCP-DCQ),0.D0)
            DB23=DCMPLX(0.D0,C15*(DSTP+DSDQ))
            DB33=DCMPLX(-C12*DCP+C11*DCQ,0.D0)
            DB43=DCMPLX(0.D0,-(C11*DSTP+C12*DSDQ))
C
            DB14=DCMPLX(0.D0,C15*(DSDP+DSTQ))
            DB24=DB13
            DB34=DCMPLX(0.D0,C12*DSDP+C11*DSTQ)
            DB44=DCMPLX(C11*DCP-C12*DCQ,0.D0)
C
C         * DCPZ - DSDQZ :
C
            DCPZ = WN*DSTP
            DCQZ = WN*DSTQ
            DSTPZ=-WN*RA2*DCP
            DSTQZ=-WN*RB2*DCQ
            DSDPZ=-WN*DCP
            DSDQZ=-WN*DCQ
C
C         * DB13Z - DB24Z *
C
            DB13Z=DCMPLX(-C15*(DCPZ-DCQZ),0.D0)
            DB23Z=DCMPLX(0.D0,C15*(DSTPZ+DSDQZ))
            DB33Z=DCMPLX(-C12*DCPZ+C11*DCQZ,0.D0)
            DB43Z=DCMPLX(0.D0,-(C11*DSTPZ+C12*DSDQZ))
C
            DB14Z=DCMPLX(0.D0,C15*(DSDPZ+DSTQZ))
            DB24Z=DB13Z
            DB34Z=DCMPLX(0.D0,C12*DSDPZ+C11*DSTQZ)
            DB44Z=DCMPLX(C11*DCPZ-C12*DCQZ,0.D0)
C
C         * L(M+1) - S(M+1) *
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)
            RG=B11*AA(1,3)+B12*AA(2,3)+B13*AA(3,3)+B14*AA(4,3)
            RR=B11*AA(1,4)+B12*AA(2,4)+B13*AA(3,4)+B14*AA(4,4)
C
            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)
            RH=B21*AA(1,3)+B22*AA(2,3)+B23*AA(3,3)+B24*AA(4,3)
            RS=B21*AA(1,4)+B22*AA(2,4)+B23*AA(3,4)+B24*AA(4,4)
C
C         * DGS - DSSZ *
C
            DGS(IL)=IMAG(RL*DB13+RK*DB23+RG*DB33+RR*DB43)
            DHS(IL)=REAL(RN*DB13+RM*DB23+RH*DB33+RS*DB43)
            DRS(IL)=REAL(RL*DB14+RK*DB24+RG*DB34+RR*DB44)
            DSS(IL)=IMAG(RN*DB14+RM*DB24+RH*DB34+RS*DB44)
C
            DGSZ(IL)=IMAG(RL*DB13Z+RK*DB23Z+RG*DB33Z+RR*DB43Z)
            DHSZ(IL)=REAL(RN*DB13Z+RM*DB23Z+RH*DB33Z+RS*DB43Z)
            DRSZ(IL)=REAL(RL*DB14Z+RK*DB24Z+RG*DB34Z+RR*DB44Z)
            DSSZ(IL)=IMAG(RN*DB14Z+RM*DB24Z+RH*DB34Z+RS*DB44Z)
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 AND ARK **
C
          UU(IOM,MODE) =-FRK/FRW
          ARK(IOM,MODE)=(RGI*RNR-RLI*RHR)/(WN*FRK)
          UW0(IOM,MODE)=-IMAG(RK/RL)
C
C      ** U0 AND W0 **
C
         HN1=0.D0
C
         DO 300 IL=1,NL-1
C
            DB1=DGS(IL)*RNR-DHS(IL)*RLI
            DB2=DRS(IL)*RNR+DSS(IL)*RLI
            DB1Z=DGSZ(IL)*RNR-DHSZ(IL)*RLI
            DB2Z=DRSZ(IL)*RNR+DSSZ(IL)*RLI
C
            DWZX=DCMPLX(DB1/FRK,0.D0)
            DWXX=DCMPLX(0.D0,-DB2/FRK)
            DWZZ=DCMPLX(0.D0,-DB1Z/(FRK*WN))
            DWXZ=DCMPLX(-DB2Z/(FRK*WN),0.D0)
C
            DUZX=-RK*DWZX/RL
            DUXX=-RK*DWXX/RL
            DUZZ=-RK*DWZZ/RL
            DUXZ=-RK*DWXZ/RL
C
            TP=2.D0/(1.D0+DNS(IL)*VP(IL)/(DS0*P0))
            TS=2.D0/(1.D0+DNS(IL)*VS(IL)/(DS0*S0))
C
            SMM=DNS(IL)*VS(IL)*VS(IL)
            RAM=DNS(IL)*VP(IL)*VP(IL)-2.D0*SMM
C
            XO=0.D0
            IF(IL.EQ.1)THEN
              ZO=0.D0
            ELSE
              ZO=ZO+DEP(IL-1)
            END IF
C
            CALL HSKSV(OM,XF,ZF,XO,ZO,P0,S0,WD,DIP,DIS,RV,RT,UF,WF)
            UF=2.D0*UF
            WF=2.D0*WF
C
            IF((IL.GE.NFSL).AND.(IL.LE.NFNL))THEN
              U0(IOM)=U0(IOM)+2.D0*TP*UF*(RAM*(DUXX+DWXZ)+2.D0*SMM*DUXX)
              U0(IOM)=U0(IOM)+2.D0*TS*WF* SMM*(DUXZ+DWXX)
              W0(IOM)=W0(IOM)+2.D0*TP*UF*(RAM*(DUZX+DWZZ)+2.D0*SMM*DUZX)
              W0(IOM)=W0(IOM)+2.D0*TS*WF* SMM*(DUZZ+DWZX)
            END IF
C
            HN1=HN1+DEP(IL)
C
  300     CONTINUE
C
          IF(NFNL.NE.NL)GOTO 110
C
          RRAN=-IMAG(RAN)
          RRBN=-IMAG(RBN)
          EXA=WN*RRAN*DDEP
C
          DUP=(0.D0,0.D0)
          DUS=(0.D0,0.D0)
          DWP=(0.D0,0.D0)
          DWS=(0.D0,0.D0)
C
          DO 400 INL=1,NNL
C
            XO=0.D0
            ZO=HN1+DDEP*DFLOAT(INL-1)
C
            CALL HSKSV(OM,XF,ZF,XO,ZO,P0,S0,WD,DIP,DIS,RV,RT,UF,WF)
            UF=2.D0*UF
            WF=2.D0*WF
C
            DP=DEXP(-EXA*DFLOAT(INL))-DEXP(-EXA*DFLOAT(INL-1))
            DS=DEXP(-EXB*DFLOAT(INL))-DEXP(-EXB*DFLOAT(INL-1))
C
            DUP=DUP+UF*DP
            DUS=DUS*UF*DS
            DWP=DWP+WF*DP
            DWS=DWS*WF*DS
C
  400     CONTINUE
C
          CCC=1.D0/(DNS(NL)*CV*CV)
          CCW=CCC/WN
          DUGN=-CCW*DUP
          DWGN=-CCW*DWP
          DUSN= CCW*DUS
          DWSN= CCW*DWS
          DUHN=-CCW*DUS/RRBN
          DWHN=-CCW*DWS/RRBN
          DURN=-CCW*DUP/RRAN
          DWRN=-CCW*DWP/RRAN
C
          DUGZ= CCC*DUP*RRAN
          DWGZ= CCC*DWP*RRAN
          DUSZ=-CCC*DUS*RRBN
          DWSZ=-CCC*DWS*RRBN
          DUHZ= CCC*DUS
          DWHZ= CCC*DWS
          DURZ= CCC*DUP
          DWRZ= CCC*DWP
C
          DUB1=DUGN*RNR-DUHN*RLI
          DWB1=DWGN*RNR-DWHN*RLI
          DUB2=DURN*RNR+DUSN*RLI
          DWB2=DWRN*RNR+DWSN*RLI
C
          DUB1Z=DUGZ*RNR-DUHZ*RLI
          DWB1Z=DWGZ*RNR-DWHZ*RLI
          DUB2Z=DURZ*RNR+DUSZ*RLI
          DWB2Z=DWRZ*RNR+DWSZ*RLI
C
          DUWZX= DUB1/FRK
          DWWZX= DWB1/FRK
          DUWXX=-DUB2/FRK*DCMPLX(0.D0,1.D0)
          DWWXX=-DWB2/FRK*DCMPLX(0.D0,1.D0)
C
          DUWZZ=-DUB1Z/(FRK*WN)*DCMPLX(0.D0,1.D0)
          DWWZZ=-DWB1Z/(FRK*WN)*DCMPLX(0.D0,1.D0)
          DUWXZ=-DUB2Z/(FRK*WN)
          DWWXZ=-DWB2Z/(FRK*WN)
C
          DUUZX=-RK*DUWZX/RL
          DWUZX=-RK*DUWZX/RL
          DUUXX=-RK*DUWXX/RL
          DWUXX=-RK*DUWXX/RL
C
          DUUZZ=-RK*DUWZZ/RL
          DWUZZ=-RK*DUWZZ/RL
          DUUXZ=-RK*DUWXZ/RL
          DWUXZ=-RK*DUWXZ/RL
C
          TP=2.D0/(1.D0+DNS(NL)*VP(NL)/(DS0*P0))
          TS=2.D0/(1.D0+DNS(NL)*VS(NL)/(DS0*S0))
C
          SMM=DNS(NL)*VS(NL)*VS(NL)
          RAM=DNS(NL)*VP(NL)*VP(NL)-2.D0*SMM
C
          U0(IOM)=U0(IOM)+2.D0*(TP*(RAM*(DUUXX+DUWXZ)+2.D0*SMM*DUUXX)
     *                         +TS* SMM*(DWUXZ+DWWXX))
          W0(IOM)=W0(IOM)+2.D0*(TP*(RAM*(DUUZX+DUWZZ)+2.D0*SMM*DUUZX)
     *                         +TS* SMM*(DWUZZ+DWWZX))
C
  110   CONTINUE
C
  100  CONTINUE
C                     
      RETURN
      END
C
C
C
      SUBROUTINE RYSC(WN,BM,CV,CB2,DEP,CQ,STQ,SDQ)
      IMPLICIT REAL*8(A-H,O-Z)
C
      IF(BM-CV)10,20,30
C
   10 RBM=DSQRT(CB2-1.D0)
      QM=WN*RBM*DEP
      CQ=DCOS(QM)
      SQ=DSIN(QM)
      STQ=RBM*SQ
      SDQ=SQ/RBM
      GOTO 40
C
   20 CQ=1.D0
      STQ=0.D0
      SDQ=WN*DEP
C
   30 RBM=DSQRT(1.D0-CB2)
      QM=WN*RBM*DEP
      EX1=DEXP(QM)
      EX2=DEXP(-QM)
      CQ=(EX1+EX2)/2.0D0
      EX12=(EX1-EX2)/2.0D0
      STQ=-RBM*EX12
      SDQ=EX12/RBM
C
   40 RETURN
      END
