C
C *** MATRIX DV AND DP DUE TO LOVE WAVE ***
C
      SUBROUTINE LVVP(OM,CC,NMODE,MMD,SZ4,SZ5,ND4,DV,DP)
C
C ** INPUT DATA **
C OM    : OMEGA
C CC(J) : PHASE VELOCITY OF Jth LOVE MODE
C NMODE : OUTPUT HIGHEST MODE
C MMD   : NEMD (HIGHEST MODE) IN MAIN PROGRAM DIMENSION
C SZ4(I): Z LOCATION OF Ith NODE ON BOUNDARY 4 (THE UPPER LAYER)
C SZ5(I): Z LOCATION OF Ith NODE ON BOUNDARY 5 (THE BOTTOM BEDROCK)
C ND4   : DIMENSION OF DV AND DP
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),SZ4(ND4),SZ5(ND4)
      COMMON /MATP/VS1,VS2,DN1,DN2,Q1,Q2 /DEPTH/DEP
     *       /NS/NS1,NS2,NS3,NS4,NS5
      COMPLEX*16 DV(ND4,ND4),DP(ND4,ND4)
     *          ,SB11,SB12,SBK11,SBK21,SBW11,SBW21,RBN
C
C *** INITIALIZATION 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 *** SHEAR MODULUS ***
C
      VS12=VS1*VS1
      VS22=VS2*VS2
      SM1=DN1*VS12
      SM2=DN2*VS22
C
C *** MODE ***
C
      DO 200 MODE=0,NMODE
C
        CV=CC(MODE)
        WN=OM/CV
C
        CV2=CV*CV
        CB2=CV2/VS12
        CALL LVSC(WN,VS1,CV,CB2,DEP,CQ,STQ,SDQ,RBM2)
C
        DSG=DEP*SDQ
        DC=DEP*CQ
        SGK=SDQ/WN
C
        SB11=DCMPLX(CQ,0.0D0)
        SB12=DCMPLX(0.0D0,SDQ/SM1)
C
        CVM2=CV/VS12
        SBK11=DCMPLX(DSG,0.0D0)
        SBW11=CVM2*DCMPLX(DSG,0.0D0)
        SBK21=DCMPLX(0.0D0,-SM1*(CB2*SGK+DC))
        SBW21=CVM2*DCMPLX(0.0D0,SM1*(SGK+DC))
C
        CB2=CV2/VS22
        RRBN=DSQRT(1.D0-CB2)
        RBN=DCMPLX(0.D0,-RRBN)
C
        FLK= -IMAG(SBK21+SM2*(RBN*SBK11-CB2*SB11/(RBN*WN)))
        FLW= -IMAG(SBW21+SM2*(RBN*SBW11+CB2*SB11/(CV*RBN*WN)))
        UU=-FLK/FLW
C
        AL11=1.D0/REAL(SM2*RBN*SB12+SB11)
        ALM=1.D0/(AL11*FLK)
C
C   *** SOURCE ON B4 ***
C
        DO 210 K=1,NS4
C
          SDEP=(SZ4(K)+SZ4(K+1))/2.D0
          CB2=CV2/VS12
          CALL LVSC(WN,VS1,CV,CB2,SDEP,CQ,STQ,SDQ,RBM2)
C
          ASM11=CQ
C
C     *** OBSERVATION ON B4 ***
C
          DO 220 J=1,NS4
C
            ODEP2=SZ4(J)
            ODEP1=SZ4(J+1)
C
            CB2=CV2/VS12
            CALL LVSC(WN,VS1,CV,CB2,ODEP2,CQ,STQ,SDQ2,RBM2)
            CALL LVSC(WN,VS1,CV,CB2,ODEP1,CQ,STQ,SDQ1,RBM2)
C
            DAOM11=(SDQ2-SDQ1)/WN
C
            DAMKJ=ALM*ASM11*DAOM11
            DVKJ=-DAMJK/WN
            DPKJ=-SM1*DAMJK
C
            DV(K,J)=DV(K,J)+DCMPLX(0.D0,DVKJ)
            DP(K,J)=DP(K,J)+DCMPLX(DPKJ,0.D0)
C
  220     CONTINUE  
C
C     *** OBSERVATION ON B5 ***
C
          DO 230 J=1,NS5
C
            ODEP2=SZ5(J+1)
            ODEP1=SZ5(J)
C
            CB2=CV2/VS22
            RRBN=DSQRT(1.D0-CB2)
            RM2=WN*RRBN
            EX2=DEXP(-RM2*(ODEP2-DEP))
            EX1=DEXP(-RM2*(ODEP1-DEP))
C
            DAOM11=-(EX2-EX1)/RM2
C
            DAMKJ=ALM*ASM11*DAOM11
            DVKJ=-DAMJK/WN
            DPKJ=-SM1*DAMJK
C
            DV(K,NS4+J)=DV(K,NS4+J)+DCMPLX(0.D0,DVKJ)
            DP(K,NS4+J)=DP(K,NS4+J)+DCMPLX(DPKJ,0.D0)
C
  230     CONTINUE  
C
  210   CONTINUE  
C
C   *** SOURCE ON B5 ***
C
        DO 250 K=1,NS5
C
          SDEP=(SZ5(K)+SZ5(K+1))/2.D0
C
          CB2=CV2/VS22
          RRBN=DSQRT(1.D0-CB2)
          RM2=WN*RRBN
          EX=DEXP(-RM2*(SDEP-DEP))
C
          ASM11=AL11*EX
C
C     *** OBSERVATION ON B4 ***
C
          DO 260 J=1,NS4
C
            ODEP2=SZ4(J)
            ODEP1=SZ4(J+1)
C
            CB2=CV2/VS12
            CALL LVSC(WN,VS1,CV,CB2,ODEP2,CQ,STQ,SDQ2,RBM2)
            CALL LVSC(WN,VS1,CV,CB2,ODEP1,CQ,STQ,SDQ1,RBM2)
C
            DAOM11=(SDQ2-SDQ1)/WN
C
            DAMKJ=ALM*ASM11*DAOM11
            DVKJ=-DAMJK/WN
            DPKJ=-SM1*DAMJK
C
            DV(NS4+K,J)=DV(NS4+K,J)+DCMPLX(0.D0,DVKJ)
            DP(NS4+K,J)=DP(NS4+K,J)+DCMPLX(DPKJ,0.D0)
C
  260     CONTINUE  
C
C     *** OBSERVATION ON B5 ***
C
          DO 270 J=1,NS5
C
            ODEP2=SZ4(J+1)
            ODEP1=SZ4(J)
C
            CB2=CV2/VS22
            RRBN=DSQRT(1.D0-CB2)
            RM2=WN*RRBN
            EX2=DEXP(-RM2*(ODEP2-DEP))
            EX1=DEXP(-RM2*(ODEP1-DEP))
C
            DAOM11=-(EX2-EX1)/RM2
C
            DAMKJ=ALM*ASM11*DAOM11
            DVKJ=-DAMJK/WN
            DPKJ=-SM1*DAMJK
C
            DV(NS4+K,NS4+J)=DV(NS4+K,NS4+J)+DCMPLX(0.D0,DVKJ)
            DP(NS4+K,NS4+J)=DP(NS4+K,NS4+J)+DCMPLX(DPKJ,0.D0)
C
  270     CONTINUE  
C
  250   CONTINUE  
C
  200 CONTINUE
C
C *** DIAGONAL TERMS ***
C
      DO 300 I=1,NS4
        DP(I,I)=DP(I,I)+(0.5D0,0.D0)
  300 CONTINUE
C        
      DO 400 I=1,NS5
        DP(NS4+I,NS4+I)=DP(NS4+I,NS4+I)+(0.5D0,0.D0)
  400 CONTINUE
C        
      RETURN
      END
C
C
C
      SUBROUTINE LVSC(WN,BM,CV,CB2,DEP,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*DEP
        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 : OM,IL,CV=',OM,IL,CV
        STOP
C
   30   RBM=DSQRT(1.0D0-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
        RBM2=-RBM*RBM
C
   40 RETURN
C
      END
