C
C *** LOVE WAVE GROUP VELOCITY AND MEDIUM RESPONSE ***
C
      SUBROUTINE LVUAL(DOM,NOM,MOM,ML,NL,DNS,VS,DEP,CC,NMODE,CFOM
     *          ,NSMD,MMD,MLS,SDEP,UU,ALM)
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    : TOTAL NUMBER OF LAYERS           
C DNS(I): DENSITY OF Ith LAYER 
C VS (I): VS OF Ith LAYER      
C DEP(I): DEPTH ( THICKNESS ) OF Ith LAYER  
C CC(I,J) : PHASE VELOCITY OF Jth MODE LOVE WAVE 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 MLS     : NUMBER OF LAYER WHERE A SOURCE IS LOCATED (FROM TOP LAYER)
C SDEP    : SOURCE DEPTH IN THE SOURCE LAYER (FROM UPPER BOUNDARY)
C ** OUTPUT DATA **
C UU(I,J) : GROUP VELOCITY OF Jth MODE LOVE WAVE AT Ith OMEGA
C ALM(I,J): AMPLITUDE SPEC (i*AL) OF Jth MODE LOVE WAVE AT Ith OMEGA
C           TO OBTAIN AN ACTUAL ALM, MULTIPLY (-i) BY AL
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION CC(MOM,0:MMD),NMODE(MOM),UU(MOM,0:MMD),ALM(MOM,0:MMD)
     *         ,CFOM(0:MMD),DNS(ML),VS(ML),DEP(ML)
      COMPLEX*16 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,AM11,AM12,AM21,AM22,BM1,BM2,FM1,FM2
C
      DO 100 IOM=1,NOM
C
        OM=DOM*DFLOAT(IOM)
        IF(OM.LT.CFOM(NSMD))GOTO 100
C
        DO 110 MODE=0,NMODE(IOM)
C
          CV=CC(IOM,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
          DO 120 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
            SM=DNS(IL)*VS(IL)*VS(IL)
            BM=VS(IL)
            CB2=CV*CV/(BM*BM)
            CALL LVSC(WN,BM,CV,CB2,DEP(IL),CQ,STQ,SDQ,RBM2)
C
            DSG=DEP(IL)*SDQ
            DC=DEP(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
            IF(IL.EQ.MLS)THEN
              AM11=SA11
              AM12=SA12
              AM21=SA21
              AM22=SA22
              CALL LVSC(WN,BM,CV,CB2,DEP(IL)-SDEP,CQ,STQ,SDQ,RBM2)
              BM1=DCMPLX(0.0D0,SDQ/SM)
              BM2=DCMPLX(CQ,0.0D0)
            END IF
C
  120     CONTINUE
C
          CB2=CV*CV/(VS(NL)*VS(NL))
          RRBN=DSQRT(1.0D0-CB2)
          RBN=DCMPLX(0.0D0,-RRBN)
C
          SM=DNS(NL)*VS(NL)*VS(NL)
          IF(MLS.EQ.NL)THEN
            SAL=DEXP(-WN*RRBN*SDEP)
          ELSE
            FM1=SM*RBN*AM11+AM21
            FM2=SM*RBN*AM12+AM22
            SAL=REAL(FM1*BM1+FM2*BM2)
          END IF
C
          FLKIM=IMAG(RK21+SM*(RBN*RK11-CB2*R11/(RBN*WN)))
          FLWIM=IMAG(RW21+SM*(RBN*RW11+CB2*R11/(CV*RBN*WN)))
          UU(IOM,MODE)=-FLKIM/FLWIM
          ALM(IOM,MODE)=SAL/FLKIM
C
  110   CONTINUE
C
  100 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
