C
C *** MATRIX DV AND DP ***
C
      SUBROUTINE DVDP(OM,ML,NL,DNS,VS,THK,CC,NMODE
     *               ,MMD,ZZ,MD4,ND4,DV,DP)
C
C ** INPUT DATA **
C OM     : OMEGA
C ML     : MAIN NL
C NL     : NUMBER OF LAYERS
C DNS(I) : DENSITY OF THE Ith LAYER 
C VS (I) : VS OF THE Ith LAYER
C THK(I) : THICKNESS OF THE Ith LAYER
C CC(J)  : PHASE VELOCITY OF THE Jth LOVE MODE
C NMODE  : OUTPUT HIGHEST MODE
C MMD    : MAX MODE NUMBER IN MAIN DIMENSION
C ZZ(I)  : DEPTH OF Ith NODE FROM THE FREE SURFACE
C MD4    : MAIN ND4
C ND4    : NUMBER OF TOTAL NODE
C ** OUTPUT DATA **
C DV(I,J): MATRIX DV
C DP(I,J): MATRIX DP
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CC(0:MMD),DNS(ML),VS(ML),THK(ML),ZZ(MD4)
      COMPLEX*16 DV(MD4,MD4),DP(MD4,MD4)
     *          ,R11,R12,R21,R22,SA11,SA12,SA21,SA22
     *          ,SB11,SB12,SB21,SB22
     *          ,RK11,RK12,RK21,RK22,SAK11,SAK12,SAK21,SAK22
     *          ,SBK11,SBK12,SBK21,SBK22
     *          ,RW11,RW12,RW21,RW22,SAW11,SAW12,SAW21,SAW22
     *          ,SBW11,SBW12,SBW21,SBW22,RBN
     *          ,R11S,R12S,R21S,R22S,R11O,R12O,R21O,R22O
     *          ,DCJ,DSJ,DZCJ,DZSJ,DVKJ,DZVKJ,DV1,DV2
C
C *** INITILIZATION OF DV AND DP ***
C
      DO 100 K=1,ND4
        DO 110 J=1,ND4
          DV(K,J)=(0.D0,0.D0)
          DP(K,J)=(0.D0,0.D0)
  110   CONTINUE
  100 CONTINUE
C
C *** LOVE MODE ***
C
      DO 200 MODE=0,NMODE
C
        CV=CC(MODE)
        WN=OM/CV
C
        R11=(1.0D0,0.0D0)
        R12=(0.0D0,0.0D0)
        R21=(0.0D0,0.0D0)
        R22=(1.0D0,0.0D0)
C
        RK11=(0.0D0,0.0D0)
        RK12=(0.0D0,0.0D0)
        RK21=(0.0D0,0.0D0)
        RK22=(0.0D0,0.0D0)
C
        RW11=(0.0D0,0.0D0)
        RW12=(0.0D0,0.0D0)
        RW21=(0.0D0,0.0D0)
        RW22=(0.0D0,0.0D0)
C
C   *** AL11 AND ALM ***
C
        DO 210 IL=NL-1, 1, -1 
C
          SA11=R11
          SA12=R12
          SA21=R21
          SA22=R22
C
          SAK11=RK11
          SAK12=RK12
          SAK21=RK21
          SAK22=RK22
C
          SAW11=RW11
          SAW12=RW12
          SAW21=RW21
          SAW22=RW22
C
          BM=VS(IL)
          SM=DNS(IL)*BM*BM
          CB2=CV*CV/(BM*BM)
          CALL LVSC(WN,BM,CV,CB2,THK(IL),CQ,STQ,SDQ,RBM2)
C
          DSG=THK(IL)*SDQ
          DC=THK(IL)*CQ
          SGK=SDQ/WN
          SG2=SM*RBM2
C
          SB11=DCMPLX(CQ,0.0D0)
          SB12=DCMPLX(0.0D0,SDQ/SM)
          SB21=DCMPLX(0.0D0,STQ*SM)
          SB22=DCMPLX(CQ,0.0D0)
C
          SBK11=DCMPLX(DSG,0.0D0)
          SBK12=DCMPLX(0.0D0,(CB2*SGK-DC)/SG2)
          SBK21=DCMPLX(0.0D0,-SM*(CB2*SGK+DC))
          SBK22=DCMPLX(DSG,0.0D0)
C
          CVM2=CV/(BM*BM)
          SBW11=CVM2*DCMPLX(-DSG,0.0D0)
          SBW12=CVM2*DCMPLX(0.0D0,(-SGK+DC)/SG2)
          SBW21=CVM2*DCMPLX(0.0D0,SM*(SGK+DC))
          SBW22=CVM2*DCMPLX(-DSG,0.0D0)
C
          R11=SA11*SB11+SA12*SB21
          R12=SA11*SB12+SA12*SB22
          R21=SA21*SB11+SA22*SB21
          R22=SA21*SB12+SA22*SB22
C
          RK11=SAK11*SB11+SAK12*SB21+SA11*SBK11+SA12*SBK21
          RK12=SAK11*SB12+SAK12*SB22+SA11*SBK12+SA12*SBK22
          RK21=SAK21*SB11+SAK22*SB21+SA21*SBK11+SA22*SBK21
          RK22=SAK21*SB12+SAK22*SB22+SA21*SBK12+SA22*SBK22
C
          RW11=SAW11*SB11+SAW12*SB21+SA11*SBW11+SA12*SBW21
          RW12=SAW11*SB12+SAW12*SB22+SA11*SBW12+SA12*SBW22
          RW21=SAW21*SB11+SAW22*SB21+SA21*SBW11+SA22*SBW21
          RW22=SAW21*SB12+SAW22*SB22+SA21*SBW12+SA22*SBW22
C
  210   CONTINUE
C
        BM=VS(NL)
        CB2=CV*CV/(BM*BM)
        RRBN=DSQRT(1.D0-CB2)
        RBN=DCMPLX(0.D0,-RRBN)
        SM=DNS(NL)*BM*BM
C
        FLK= -IMAG(RK21+SM*(RBN*RK11-CB2*R11/(RBN*WN)))
        FLW= -IMAG(RW21+SM*(RBN*RW11+CB2*R11/(CV*RBN*WN)))
        UU=-FLK/FLW
C
        AL11=1.D0/REAL(SM*RBN*R12+R22)
        ALM=1.D0/(AL11*FLK)
C
C   *** DEPTH OF TOP OF BOTTOM BEDROCK FROM FREE SURFACE ***
C
        DNL=0.D0
        DO 220 K=1,NL-1
          DNL=DNL+THK(K)
  220   CONTINUE
C
C   *** SOURCE POINT ***
C
        R11S=(1.0D0,0.0D0)
        R12S=(0.0D0,0.0D0)
        R21S=(0.0D0,0.0D0)
        R22S=(1.0D0,0.0D0)
C
        ILS=1
        HDEP1S=0.D0
        HDEP2S=THK(ILS)
C
        DO 230 NSP=1,ND4
C
          SDEP=ZZ(NSP)
          IF(ILS.EQ.NL)GOTO 10
C
          IF(SDEP.GE.HDEP2S)THEN
C
            SA11=R11S
            SA12=R12S
            SA21=R21S
            SA22=R22S
C
            BM=VS(ILS)
            SM=DNS(ILS)*BM*BM
            CB2=CV*CV/(BM*BM)
            CALL LVSC(WN,BM,CV,CB2,THK(ILS),CQ,STQ,SDQ,RBM2)
C
            SB11=DCMPLX(CQ,0.0D0)
            SB12=DCMPLX(0.0D0,SDQ/SM)
            SB21=DCMPLX(0.0D0,STQ*SM)
            SB22=DCMPLX(CQ,0.0D0)
C
            R11S=SB11*SA11+SB12*SA21
            R12S=SB11*SA12+SB12*SA22
            R21S=SB21*SA11+SB22*SA21
            R22S=SB21*SA12+SB22*SA22
C
            ILS=ILS+1
            HDEP1S=HDEP2S
            IF(ILS.EQ.NL)THEN
              HDEP2S=HDEP2S*1000.D0
            ELSE
              HDEP2S=HDEP2S+THK(ILS)
            END IF
C
          END IF
C
   10     IF(ILS.EQ.NL)THEN
C
            BM=VS(NL)
            CB2=CV*CV/(BM*BM)
            RRBN=DSQRT(1.D0-CB2)
            RNL=WN*RRBN
C
            ASM11=AL11*DEXP(-RNL*(SDEP-DNL))
C
          ELSE
C
            BM=VS(ILS)
            SM=DNS(ILS)*BM*BM
            CB2=CV*CV/(BM*BM)
            CALL LVSC(WN,BM,CV,CB2,SDEP-HDEP1S,CQ,STQ,SDQ,RBM2)
C
            SBK11=DCMPLX(CQ,0.0D0)
            SBK12=DCMPLX(0.0D0,STQ/SM)
C
            ASM11=REAL(SBK11*R11S+SBK12*R21S)
C
          END IF
C
C     *** OBSERVATION POINT ***
C
          R11O=(1.0D0,0.0D0)
          R12O=(0.0D0,0.0D0)
          R21O=(0.0D0,0.0D0)
          R22O=(1.0D0,0.0D0)
C
          ILO=1
          HDEP1O=0.D0
          HDEP2O=THK(ILO)
C
          DO 240 NOP=1,ND4-1
C
            ODEP1=ZZ(NOP)
            ODEP2=ZZ(NOP+1)
            IF(ILO.EQ.NL)GOTO 20
C
            IF(ODEP1.GE.HDEP2O)THEN
C
              SA11=R11O
              SA12=R12O
              SA21=R21O
              SA22=R22O
C
              BM=VS(ILO)
              SM=DNS(ILO)*BM*BM
              CB2=CV*CV/(BM*BM)
              CALL LVSC(WN,BM,CV,CB2,THK(ILO),CQ,STQ,SDQ,RBM2)
C
              SB11=DCMPLX(CQ,0.0D0)
              SB12=DCMPLX(0.0D0,SDQ/SM)
              SB21=DCMPLX(0.0D0,STQ*SM)
              SB22=DCMPLX(CQ,0.0D0)
C
              R11O=SB11*SA11+SB12*SA21
              R12O=SB11*SA12+SB12*SA22
              R21O=SB21*SA11+SB22*SA21
              R22O=SB21*SA12+SB22*SA22
C
              ILO=ILO+1
              HDEP1O=HDEP2O
              IF(ILO.EQ.NL)THEN
                HDEP2O=HDEP2O*1000.D0
              ELSE
                HDEP2O=HDEP2O+THK(ILO)
              END IF
C
            END IF
C
   20       IF(ILO.EQ.NL)THEN
C
              BM=VS(NL)
              CB2=CV*CV/(BM*BM)
              RRBN=DSQRT(1.D0-CB2)
              RNL=WN*RRBN
C
              EX1=DEXP(-RNL*(ODEP1-DNL))
              EX2=DEXP(-RNL*(ODEP2-DNL))
              DEM=-(EX2-EX1)/RNL
              DZEM=(-(RNL*SDEP2+1)*EX2+(RNL*SDEP1+1)*EX1)/(RNL*RNL)
C
              DAOM11=AL11*DEM
              DZAOM11=AL11*DZEM
C
            ELSE
C
              BM=VS(ILO)
              SM=DNS(ILO)*BM*BM
              CB2=CV*CV/(BM*BM)
              CALL LVSC(WN,BM,CV,CB2,ODEP2-HDEP1O
     *                 ,CQ2,STQ2,SDQ2,RBM2)
              CALL LVSC(WN,BM,CV,CB2,ODEP1-HDEP1O
     *                 ,CQ1,STQ1,SDQ1,RBM2)
C
              DCJ=DCMPLX(SDQ2-SDQ1,0.0D0)/WN
              DSJ=DCMPLX(0.D0,-(CQ2-CQ1)/(WN*SM*RBM2))
C
              DZC2=CQ2+WN*SDEP2*STQ2
              DZC1=CQ1+WN*SDEP1*STQ1
              DZCJ=DCMPLX((DZC2-DZC1)/(WN*WN*RBM2),0.0D0)
              DZS2=SDQ2-WN*SDEP2*CQ2
              DZS1=SDQ1-WN*SDEP1*CQ1 
              DZSJ=DCMPLX(0.D0,(DZS2-DZS1)/(SM*WN*WN*RBM2))
C
              DAOM11=REAL(DCJ*R11O+DSJ*R21)
              DZAOM11=REAL(DZCJ*R11O+DZSJ*R21)
C
            END IF
C
            DAKJ=ALM*ASM11*DAOM11
            DZAKJ=ALM*ASM11*DZAOM11
C
            SM=DNS(ILO)*VS(ILO)*VS(ILO)
            DVKJ=-DCMPLX(0.D0,1.D0)*DAKJ/WN
            DPKJ=-SM*DAKJ
            DZVKJ=-DCMPLX(0.D0,1.D0)*DZAKJ/WN
            DZPKJ=-SM*DZAKJ
C
            DH=1.D0+ODEP1/THK(ILO)
            DV1=DVKJ*DH-DZVJK
            DV2=(DZVJK-DVJK*ODEP1)/THK(ILO)
            DP1=DPKJ*DH-DZPJK
            DP2=(DZPJK-DPJK*ODEP1)/THK(ILO)
C
            DV(NSP,NOP)=DV(NSP,NOP)+DV1
            DV(NSP,NOP+1)=DV(NSP,NOP+1)+DV2
            DP(NSP,NOP)=DP(NSP,NOP)+DCMPLX(DP1,0.D0)
            DP(NSP,NOP+1)=DP(NSP,NOP+1)+DCMPLX(DP2,0.D0)
C
  240     CONTINUE
C
  230   CONTINUE
C
  200 CONTINUE
C
C *** DIAGONAL TERM OF DP ***
C
      DO 300 I=1,ND4
        DP(I,I)=DP(I,I)+(0.5D0,0.D0)
  300 CONTINUE
C
      WRITE(*,*)
C
      RETURN
      END
C
C
C
      SUBROUTINE LVSC(WN,BM,CV,CB2,THK,CQ,STQ,SDQ,RBM2)
      IMPLICIT REAL*8(A-H,O-Z)
C
      IF(BM-CV)10,20,30
C
   10   RBM=DSQRT(CB2-1.0D0)
        QM=WN*RBM*THK
        CQ=DCOS(QM)
        SQ=DSIN(QM)
        STQ=RBM*SQ
        SDQ=SQ/RBM
        RBM2=RBM*RBM
        GOTO 40
C
   20   WRITE(*,*)' FROM SUBROUTINE LVSC, LVUAM :'
        WRITE(*,*)' RBM=0 THEN STOP !!!!'
        WRITE(*,*)' DATA : WN'WN
        STOP
C
   30   RBM=DSQRT(1.0D0-CB2)
        QM=WN*RBM*THK
        EX1=DEXP(QM)
        EX2=DEXP(-QM)
        CQ=(EX1+EX2)/2.0D0
        EX12=(EX1-EX2)/2.0D0
        STQ=-RBM*EX12
        SDQ=EX12/RBM
        RBM2=-RBM*RBM
C
   40 RETURN
C
      END
